• Main Page
  • Namespaces
  • Classes
  • Files

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

Go to the documentation of this file.
00001 
00002 
00003 #include <CompuCell3D/Simulator.h>
00004 #include <CompuCell3D/Automaton/Automaton.h>
00005 #include <CompuCell3D/Potts3D/Potts3D.h>
00006 #include <CompuCell3D/Potts3D/CellInventory.h>
00007 #include <CompuCell3D/Field3D/WatchableField3D.h>
00008 #include <CompuCell3D/Field3D/Field3DImpl.h>
00009 #include <CompuCell3D/Field3D/Field3D.h>
00010 #include <CompuCell3D/Field3D/Field3DIO.h>
00011 #include <BasicUtils/BasicClassGroup.h>
00012 #include <CompuCell3D/steppables/BoxWatcher/BoxWatcher.h>
00013 
00014 
00015 #include <BasicUtils/BasicString.h>
00016 #include <BasicUtils/BasicException.h>
00017 #include <BasicUtils/BasicRandomNumberGenerator.h>
00018 #include <PublicUtilities/StringUtils.h>
00019 #include <string>
00020 #include <cmath>
00021 #include <iostream>
00022 #include <fstream>
00023 #include <sstream>
00024 
00025 #include <time.h>
00026 
00027 #define EXP_STL
00028 #include "FastDiffusionSolver2DFE.h"
00029 
00031 // std::ostream & operator<<(std::ostream & out,CompuCell3D::DiffusionData & diffData){
00032 //    
00033 //    
00034 // }
00035 
00036 
00037 using namespace CompuCell3D;
00038 using namespace std;
00039 
00040 
00042 void FastDiffusionSolver2DSerializer::serialize(){
00043 
00044         for(int i = 0 ; i < solverPtr->diffSecrFieldTuppleVec.size() ; ++i){
00045                 ostringstream outName;
00046 
00047                 outName<<solverPtr->diffSecrFieldTuppleVec[i].diffData.fieldName<<"_"<<currentStep<<"."<<serializedFileExtension;
00048                 ofstream outStream(outName.str().c_str());
00049                 solverPtr->outputField( outStream,solverPtr->concentrationFieldVector[i]);
00050         }
00051 
00052 }
00053 
00054 
00056 void FastDiffusionSolver2DSerializer::readFromFile(){
00057         try{
00058                 for(int i = 0 ; i < solverPtr->diffSecrFieldTuppleVec.size() ; ++i){
00059                         ostringstream inName;
00060                         inName<<solverPtr->diffSecrFieldTuppleVec[i].diffData.fieldName<<"."<<serializedFileExtension;
00061 
00062                         solverPtr->readConcentrationField(inName.str().c_str(),solverPtr->concentrationFieldVector[i]);;
00063                 }
00064 
00065         } catch (BasicException &e) {
00066                 cerr<<"COULD NOT FIND ONE OF THE FILES"<<endl;
00067                 throw BasicException("Error in reading diffusion fields from file",e);
00068         }
00069 
00070 
00071 }
00072 
00073 
00074 
00076 FastDiffusionSolver2DFE::FastDiffusionSolver2DFE()
00077 : DiffusableVector2D<float>(),deltaX(1.0),deltaT(1.0)
00078 {
00079         serializerPtr=0;
00080         serializeFlag=false;
00081         readFromFileFlag=false;
00082         haveCouplingTerms=false;
00083         serializeFrequency=0;
00084         boxWatcherSteppable=0;
00085         //    useBoxWatcher=false;
00086 
00087 }
00089 FastDiffusionSolver2DFE::~FastDiffusionSolver2DFE()
00090 {
00091 
00092         if(serializerPtr)
00093                 delete serializerPtr ; serializerPtr=0;
00094 }
00096 void FastDiffusionSolver2DFE::init(Simulator *simulator, CC3DXMLElement *_xmlData) {
00097 
00098 
00099 
00100 
00101         simPtr=simulator;
00102         potts = simulator->getPotts();
00103         automaton=potts->getAutomaton();
00104 
00106         cellInventoryPtr=& potts->getCellInventory(); 
00107 
00109         cellFieldG=(WatchableField3D<CellG *> *)potts->getCellFieldG();
00110         fieldDim=cellFieldG->getDim();
00111 
00112         update(_xmlData,true);
00113 
00114 
00115 
00117         diffusePtr=&FastDiffusionSolver2DFE::diffuse;
00118         secretePtr=&FastDiffusionSolver2DFE::secrete;
00119 
00120 
00121         numberOfFields=diffSecrFieldTuppleVec.size();
00122 
00123 
00124         vector<string> concentrationFieldNameVectorTmp; //temporary vector for field names
00126         concentrationFieldNameVectorTmp.assign(diffSecrFieldTuppleVec.size(),string(""));
00127 
00128         cerr<<"diffSecrFieldTuppleVec.size()="<<diffSecrFieldTuppleVec.size()<<endl;
00129 
00130         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00131                 //       cerr<<" concentrationFieldNameVector[i]="<<diffDataVec[i].fieldName<<endl;
00132                 //       concentrationFieldNameVector.push_back(diffDataVec[i].fieldName);
00133                 concentrationFieldNameVectorTmp[i] = diffSecrFieldTuppleVec[i].diffData.fieldName;
00134                 cerr<<" concentrationFieldNameVector[i]="<<concentrationFieldNameVectorTmp[i]<<endl;
00135         }
00136 
00137         //setting up couplingData - field-field interaction terms
00138         vector<CouplingData>::iterator pos;
00139 
00140         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00141                 pos=diffSecrFieldTuppleVec[i].diffData.couplingDataVec.begin();
00142                 for(int j = 0 ; j < diffSecrFieldTuppleVec[i].diffData.couplingDataVec.size() ; ++j){
00143 
00144                         for(int idx=0; idx<concentrationFieldNameVectorTmp.size() ; ++idx){
00145                                 if( concentrationFieldNameVectorTmp[idx] == diffSecrFieldTuppleVec[i].diffData.couplingDataVec[j].intrFieldName ){
00146                                         diffSecrFieldTuppleVec[i].diffData.couplingDataVec[j].fieldIdx=idx;
00147                                         haveCouplingTerms=true; //if this is called at list once we have already coupling terms and need to proceed differently with scratch field initialization
00148                                         break;
00149                                 }
00150                                 //this means that required interacting field name has not been found
00151                                 if( idx == concentrationFieldNameVectorTmp.size()-1 ){
00152                                         //remove this interacting term
00153                                         //                pos=&(diffDataVec[i].degradationDataVec[j]);
00154                                         diffSecrFieldTuppleVec[i].diffData.couplingDataVec.erase(pos);
00155                                 }
00156                         }
00157                         ++pos;
00158                 }
00159         } 
00160 
00161 
00162         cerr<<"FIELDS THAT I HAVE"<<endl; 
00163         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00164                 cerr<<"Field "<<i<<" name: "<<concentrationFieldNameVectorTmp[i]<<endl;
00165         }
00166 
00167 
00168 
00169 
00170         cerr<<"FlexibleDiffusionSolverFE: extra Init in read XML"<<endl;
00171 
00172         workFieldDim=Dim3D(fieldDim.x+2,fieldDim.y+2,fieldDim.z+2);
00174         if(!haveCouplingTerms){
00175                 allocateDiffusableFieldVector(diffSecrFieldTuppleVec.size()+1,workFieldDim); //+1 is for additional scratch field
00176         }else{
00177                 allocateDiffusableFieldVector(2*diffSecrFieldTuppleVec.size(),workFieldDim); //with coupling terms every field need to have its own scratch field
00178         }
00179 
00180         //here I need to copy field names from concentrationFieldNameVectorTmp to concentrationFieldNameVector
00181         //because concentrationFieldNameVector is reallocated with default values once I call allocateDiffusableFieldVector
00182 
00183 
00184         for(unsigned int i=0 ; i < concentrationFieldNameVectorTmp.size() ; ++i){
00185                 concentrationFieldNameVector[i]=concentrationFieldNameVectorTmp[i];
00186         }
00187 
00188 
00189         //register fields once they have been allocated
00190         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00191                 simPtr->registerConcentrationField(concentrationFieldNameVector[i] , concentrationFieldVector[i]);
00192                 cerr<<"registring field: "<<concentrationFieldNameVector[i]<<" field address="<<concentrationFieldVector[i]<<endl;
00193         }
00194 
00195 
00196 
00197 
00198 
00199         //    exit(0);
00200 
00201         periodicBoundaryCheckVector.assign(3,false);
00202         string boundaryName;
00203         boundaryName=potts->getBoundaryXName();
00204         changeToLower(boundaryName);
00205         if(boundaryName=="periodic")  {
00206                 periodicBoundaryCheckVector[0]=true; 
00207         }
00208         boundaryName=potts->getBoundaryYName();
00209         changeToLower(boundaryName);
00210         if(boundaryName=="periodic")  {
00211                 periodicBoundaryCheckVector[1]=true;
00212         }
00213 
00214         boundaryName=potts->getBoundaryZName();
00215         changeToLower(boundaryName);
00216         if(boundaryName=="periodic")  {
00217                 periodicBoundaryCheckVector[2]=true; 
00218         }
00219 
00220 
00221         simulator->registerSteerableObject(this);
00222 
00223 }
00225 void FastDiffusionSolver2DFE::extraInit(Simulator *simulator){
00226 
00227         if((serializeFlag || readFromFileFlag) && !serializerPtr){
00228                 serializerPtr=new FastDiffusionSolver2DSerializer();
00229                 serializerPtr->solverPtr=this;
00230         }
00231 
00232         if(serializeFlag){
00233                 simulator->registerSerializer(serializerPtr);
00234         }
00235 
00236         //checking if box watcher is necessary at all
00237 
00238 
00239         bool useBoxWatcher=false;
00240         for (int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00241                 if(diffSecrFieldTuppleVec[i].diffData.useBoxWatcher){
00242                         useBoxWatcher=true;
00243                         break;
00244                 }
00245         }
00246         bool steppableAlreadyRegisteredFlag;
00247         if(useBoxWatcher){
00248                 boxWatcherSteppable=(BoxWatcher*)Simulator::steppableManager.get("BoxWatcher",&steppableAlreadyRegisteredFlag);
00249                 if(!steppableAlreadyRegisteredFlag)
00250                         boxWatcherSteppable->init(simulator);
00251         }
00252 }
00253 
00255 void FastDiffusionSolver2DFE::start() {
00256         //     if(diffConst> (1.0/6.0-0.05) ){ //hard coded condtion for stability of the solutions - assume dt=1 dx=dy=dz=1
00257         // 
00258         //       cerr<<"CANNOT SOLVE DIFFUSION EQUATION: STABILITY PROBLEM - DIFFUSION CONSTANT TOO LARGE. EXITING..."<<endl;
00259         //       exit(0);
00260         // 
00261         //    }
00262 
00263 
00264         dt_dx2=deltaT/(deltaX*deltaX);
00265         if(readFromFileFlag){
00266                 try{
00267 
00268                         //          serializerPtr->readFromFile();
00269 
00270                 } catch (BasicException &e){
00271                         cerr<<"Going to fail-safe initialization"<<endl;
00272                         initializeConcentration(); //if there was error, initialize using failsafe defaults
00273                 }
00274 
00275         }else{
00276                 initializeConcentration();//Normal reading from User specified files
00277         }
00278 
00279 }
00280 
00282 
00283 void FastDiffusionSolver2DFE::initializeConcentration(){
00284 
00285 
00286         for(unsigned int i = 0 ; i <diffSecrFieldTuppleVec.size() ; ++i){
00287                 if(diffSecrFieldTuppleVec[i].diffData.concentrationFileName.empty()) continue;
00288                 cerr<<"fail-safe initialization "<<diffSecrFieldTuppleVec[i].diffData.concentrationFileName<<endl;
00289                 readConcentrationField(diffSecrFieldTuppleVec[i].diffData.concentrationFileName,concentrationFieldVector[i]);
00290         }
00291 
00292 
00293 }
00294 
00296 void FastDiffusionSolver2DFE::step(const unsigned int _currentStep) {
00297 
00298         currentStep=_currentStep;
00299 
00300         (this->*secretePtr)();
00301 
00302         (this->*diffusePtr)();
00303 
00304 
00305         if(serializeFrequency>0 && serializeFlag && !(_currentStep % serializeFrequency)){
00306                 serializerPtr->setCurrentStep(currentStep);
00307                 serializerPtr->serialize();
00308         }
00309 
00310 }
00312 
00313 void FastDiffusionSolver2DFE::secreteOnContactSingleField(unsigned int idx){
00314 
00315 
00316         SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;
00317 
00318 
00319         std::map<unsigned char,SecretionOnContactData>::iterator mitr;
00320         std::map<unsigned char,SecretionOnContactData>::iterator end_mitr=secrData.typeIdSecrOnContactDataMap.end();
00321 
00322         CellG *currentCellPtr;
00323 
00324         ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
00325         Array2D_t & concentrationArray = concentrationFieldPtr->getContainer();
00326 
00327         float currentConcentration;
00328         float secrConst;
00329         float secrConstMedium=0.0;
00330         std::map<unsigned char, float> * contactCellMapMediumPtr;
00331         std::map<unsigned char, float> * contactCellMapPtr;
00332         std::map<unsigned char, float>::iterator mitrTypeConst;
00333 
00334         bool secreteInMedium=false;
00335         //the assumption is that medium has type ID 0
00336         if(secrData.secretionOnContactTypeIds.find(automaton->getTypeId("Medium")) != secrData.secretionOnContactTypeIds.end()){ 
00337                 mitr=secrData.typeIdSecrOnContactDataMap.find(automaton->getTypeId("Medium"));
00338                 secreteInMedium=true;
00339                 contactCellMapMediumPtr = &(mitr->second.contactCellMap);
00340         }
00341 
00342 
00343 
00344         Point3D pt;
00345         Neighbor n;
00346         CellG *nCell=0;
00347         WatchableField3D<CellG *> *fieldG =(WatchableField3D<CellG *> *) potts->getCellFieldG();
00348         unsigned char type;
00349 
00350 
00351 
00352         for (int x = 1; x < workFieldDim.x-1; x++)
00353                 for (int y = 1; y < workFieldDim.y-1; y++){
00354                         pt=Point3D(x-1,y-1,0);
00355                         currentCellPtr=cellFieldG->get(pt);
00356                         currentConcentration = concentrationArray[x][y];
00357 
00358 
00359                         if(secreteInMedium && ! currentCellPtr){
00360 
00361                                 for (int i = 0  ; i<offsetVec.size() ; ++i ){ 
00362                                         n=boundaryStrategy->getNeighborDirect(pt,i);
00363                                         if(!n.distance)//not a valid neighbor
00364                                                 continue;
00365                                         nCell = fieldG->get(n.pt);
00366 
00367                                         if(nCell)
00368                                                 type=nCell->type;
00369                                         else
00370                                                 type=0;
00371 
00372 
00373                                         mitrTypeConst=contactCellMapMediumPtr->find(type);
00374 
00375                                         if(mitrTypeConst != contactCellMapMediumPtr->end()){//OK to secrete, contact detected
00376                                                 secrConstMedium = mitrTypeConst->second;
00377 
00378                                                 concentrationArray[x][y]=currentConcentration+secrConstMedium;
00379                                         }                     
00380 
00381 
00382 
00383                                 }                  
00384 
00385                                 continue;
00386                         }
00387 
00388                         if(currentCellPtr){
00389                                 if(secrData.secretionOnContactTypeIds.find(currentCellPtr->type) == secrData.secretionOnContactTypeIds.end()){  
00390                                     continue;
00391                                 }
00392 
00393                                 mitr=secrData.typeIdSecrOnContactDataMap.find(currentCellPtr->type);
00394                                 if(mitr!=end_mitr){
00395 
00396                                         contactCellMapPtr = &(mitr->second.contactCellMap);
00397 
00398                                         for (int i = 0  ; i<offsetVec.size() ; ++i ){ 
00399 
00400                                                 n=boundaryStrategy->getNeighborDirect(pt,i);
00401                                                 if(!n.distance)//not a valid neighbor
00402                                                         continue;
00403 
00404                                                 nCell = fieldG->get(n.pt);
00405 
00406                                                 if(nCell)
00407                                                         type=nCell->type;
00408                                                 else
00409                                                         type=0;
00410 
00411                                                 mitrTypeConst=contactCellMapPtr->find(type);
00412                                                 if(mitrTypeConst != contactCellMapPtr->end()){//OK to secrete, contact detected
00413                                                         secrConst=mitrTypeConst->second;
00414                                                         //                         concentrationField->set(pt,currentConcentration+secrConst);
00415                                                         concentrationArray[x][y]=currentConcentration+secrConst;
00416                                                 }
00417 
00418                                         }
00419 
00420 
00421                                 }
00422                         }
00423                 }
00424 
00425 }
00426 
00428 void FastDiffusionSolver2DFE::secreteSingleField(unsigned int idx){
00429 
00430 
00431         SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;
00432 
00433 
00434         std::map<unsigned char,float>::iterator mitr;   
00435         std::map<unsigned char,float>::iterator end_mitr=secrData.typeIdSecrConstMap.end();
00436         std::map<unsigned char,UptakeData>::iterator mitrUptake;
00437         std::map<unsigned char,UptakeData>::iterator end_mitrUptake=secrData.typeIdUptakeDataMap.end();
00438 
00439 
00440 
00441         CellG *currentCellPtr;
00442         Field3DImpl<float> * concentrationField=concentrationFieldVector[idx];
00443         float currentConcentration;
00444         float secrConst;
00445         float secrConstMedium=0.0;
00446         float maxUptakeInMedium=0.0;
00447         float relativeUptakeRateInMedium=0.0;
00448 
00449         ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
00450         Array2D_t & concentrationArray = concentrationFieldPtr->getContainer();
00451 
00452 
00453         bool doUptakeFlag=false;
00454         bool uptakeInMediumFlag=false;
00455         bool secreteInMedium=false;
00456         //the assumption is that medium has type ID 0
00457         mitr=secrData.typeIdSecrConstMap.find(automaton->getTypeId("Medium"));
00458 
00459         if( mitr != end_mitr){
00460                 secreteInMedium=true;
00461                 secrConstMedium=mitr->second;
00462         }
00463 
00464         //uptake for medium setup
00465         if(secrData.typeIdUptakeDataMap.size()){
00466                 doUptakeFlag=true;
00467         }
00468         //uptake for medium setup
00469         if(doUptakeFlag){
00470                 mitrUptake=secrData.typeIdUptakeDataMap.find(automaton->getTypeId("Medium"));
00471                 if(mitrUptake != end_mitrUptake){
00472                         maxUptakeInMedium=mitrUptake->second.maxUptake;
00473                         relativeUptakeRateInMedium=mitrUptake->second.relativeUptakeRate;
00474                         uptakeInMediumFlag=true;
00475 
00476                 }
00477         }
00478 
00479         Point3D pt;
00480         for (int x = 1; x < workFieldDim.x-1; x++)
00481                 for (int y = 1; y < workFieldDim.y-1; y++){
00482                         pt=Point3D(x-1,y-1,0);
00483                         currentCellPtr=cellFieldG->get(pt);
00484                         currentConcentration = concentrationArray[x][y];
00485 
00486                         if(secreteInMedium && ! currentCellPtr){
00487                                 concentrationArray[x][y]=currentConcentration+secrConstMedium;
00488                         }
00489 
00490                         if(currentCellPtr){
00491                                  if(secrData.secretionTypeIds.find(currentCellPtr->type) == secrData.secretionTypeIds.end()){ 
00492                                     continue;
00493                                  }
00494 
00495                                 mitr=secrData.typeIdSecrConstMap.find(currentCellPtr->type);
00496                                 if(mitr!=end_mitr){
00497                                         secrConst=mitr->second;
00498                                         concentrationArray[x][y]=currentConcentration+secrConst;
00499                                 }
00500                         }
00501                                 if(doUptakeFlag){
00502                                         if(uptakeInMediumFlag && ! currentCellPtr){
00503                                                 if(currentConcentration>maxUptakeInMedium){
00504                                                         concentrationArray[x][y]-=maxUptakeInMedium;    
00505                                                 }else{
00506                                                         concentrationArray[x][y]-= currentConcentration*relativeUptakeRateInMedium;     
00507                                                 }
00508                                         }
00509                                         if(currentCellPtr){
00510                                                 
00511                                                 mitrUptake=secrData.typeIdUptakeDataMap.find(currentCellPtr->type);
00512                                                 if(mitrUptake!=end_mitrUptake){
00513                                                         if(currentConcentration > mitrUptake->second.maxUptake){
00514                                                                 concentrationArray[x][y]-=mitrUptake->second.maxUptake;
00515                                                                 //cerr<<" uptake concentration="<< currentConcentration<<" relativeUptakeRate="<<mitrUptake->second.relativeUptakeRate<<" subtract="<<mitrUptake->second.maxUptake<<endl;
00516                                                         }else{
00517                                                                 //cerr<<"concentration="<< currentConcentration<<" relativeUptakeRate="<<mitrUptake->second.relativeUptakeRate<<" subtract="<<currentConcentration*mitrUptake->second.relativeUptakeRate<<endl;
00518                                                                 concentrationArray[x][y]-= currentConcentration*mitrUptake->second.relativeUptakeRate;  
00519                                                         }
00520                                                 }
00521                                         }
00522                                 }
00523                 }
00524 
00525 
00526 }
00527 
00529 void FastDiffusionSolver2DFE::secreteConstantConcentrationSingleField(unsigned int idx){
00530 
00531 
00532         SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;
00533 
00534 
00535         std::map<unsigned char,float>::iterator mitr;
00536         std::map<unsigned char,float>::iterator end_mitr=secrData.typeIdSecrConstMap.end();
00537 
00538         CellG *currentCellPtr;
00539         Field3DImpl<float> * concentrationField=concentrationFieldVector[idx];
00540         float currentConcentration;
00541         float secrConst;
00542         float secrConstMedium=0.0;
00543 
00544         ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
00545         Array2D_t & concentrationArray = concentrationFieldPtr->getContainer();
00546 
00547 
00548         bool secreteInMedium=false;
00549         //the assumption is that medium has type ID 0
00550         if(secrData.constantConcentrationTypeIds.find(automaton->getTypeId("Medium")) != secrData.constantConcentrationTypeIds.end()){
00551                 mitr=secrData.typeIdSecrConstMap.find(automaton->getTypeId("Medium"));
00552                 secreteInMedium=true;
00553                 secrConstMedium=mitr->second;
00554             
00555         }
00556 
00557 
00558         Point3D pt;
00559         for (int x = 1; x < workFieldDim.x-1; x++)
00560                 for (int y = 1; y < workFieldDim.y-1; y++){
00561                         pt=Point3D(x-1,y-1,0);
00562                         currentCellPtr=cellFieldG->get(pt);
00563                         //currentConcentration = concentrationArray[x][y];
00564 
00565                         if(secreteInMedium && ! currentCellPtr){
00566                                 concentrationArray[x][y]=secrConstMedium;
00567                         }
00568 
00569                         if(currentCellPtr){
00570                                  if(secrData.constantConcentrationTypeIds.find(currentCellPtr->type) == secrData.constantConcentrationTypeIds.end()){ 
00571                                  continue;
00572                                  }
00573 
00574                                 mitr=secrData.typeIdSecrConstMap.find(currentCellPtr->type);
00575                                 if(mitr!=end_mitr){
00576                                         secrConst=mitr->second;
00577                                         concentrationArray[x][y]=secrConst;
00578                                 }
00579                         }
00580                 }
00581 
00582 
00583 }
00584 
00585 
00586 
00588 void FastDiffusionSolver2DFE::secrete() {
00589 
00590         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){
00591 
00592                 for(unsigned int j = 0 ; j <diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec.size() ; ++j){
00593                         (this->*diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec[j])(i);
00594                 }
00595 
00596         }
00597 
00598 
00599 
00600 }
00601 
00603 float FastDiffusionSolver2DFE::couplingTerm(Point3D & _pt,std::vector<CouplingData> & _couplDataVec,float _currentConcentration){
00604 
00605         float couplingTerm=0.0;
00606         float coupledConcentration;
00607         for(int i =  0 ; i < _couplDataVec.size() ; ++i){
00608                 coupledConcentration=concentrationFieldVector[_couplDataVec[i].fieldIdx]->get(_pt);
00609                 couplingTerm+=_couplDataVec[i].couplingCoef*_currentConcentration*coupledConcentration;
00610         }
00611 
00612         return couplingTerm;
00613 
00614 
00615 
00616 }
00618 void FastDiffusionSolver2DFE::boundaryConditionInit(ConcentrationField_t *concentrationField){
00619         Array2D_t & concentrationArray=concentrationField->getContainer();
00620 
00621 
00622         if(periodicBoundaryCheckVector[0]){//periodic BC along X
00623                 for(int y=0 ; y<workFieldDim.y ; ++y){
00624                         concentrationArray[0][y]=concentrationArray[workFieldDim.x-2][y];
00625                         concentrationArray[workFieldDim.x-1][y]=concentrationArray[1][y];
00626                 }
00627 
00628 
00629         }else{//noFlux BC along X
00630                 for(int y=0 ; y<workFieldDim.y ; ++y){
00631                         concentrationArray[0][y]=concentrationArray[1][y];
00632                         concentrationArray[workFieldDim.x-1][y]=concentrationArray[workFieldDim.x-2][y];
00633                 }
00634 
00635 
00636         }
00637 
00638         if(periodicBoundaryCheckVector[1]){//periodic BC along Y
00639                 for(int x=0 ; x<workFieldDim.x ; ++x){
00640                         concentrationArray[x][0]=concentrationArray[x][workFieldDim.y-2];
00641                         concentrationArray[x][workFieldDim.y-1]=concentrationArray[x][1];
00642                 }
00643 
00644         }else{//noFlux BC along Y
00645                 for(int x=0 ; x<workFieldDim.x ; ++x){
00646                         concentrationArray[x][0]=concentrationArray[x][1];
00647                         concentrationArray[x][workFieldDim.y-1]=concentrationArray[x][workFieldDim.y-2];
00648                 }
00649 
00650         }
00651 
00652         //    if(periodicBoundaryCheckVector[0] || periodicBoundaryCheckVector[1]){
00653         //       for(int y=0 ; y<workFieldDim.y ; ++y){
00654         //             concentrationArray[0][y]=concentrationArray[workFieldDim.x-2][y];
00655         //             concentrationArray[workFieldDim.x-1][y]=concentrationArray[1][y];
00656         //       }
00657         // 
00658         //       for(int x=0 ; x<workFieldDim.x ; ++x){
00659         //             concentrationArray[x][0]=concentrationArray[x][workFieldDim.y-2];
00660         //             concentrationArray[x][workFieldDim.y-1]=concentrationArray[x][1];
00661         //       }
00662         // 
00663         //       
00664         //    }else{
00665         //       for(int y=0 ; y<workFieldDim.y ; ++y){
00666         //             concentrationArray[0][y]=concentrationArray[1][y];
00667         //             concentrationArray[workFieldDim.x-1][y]=concentrationArray[workFieldDim.x-2][y];
00668         //       }
00669         // 
00670         //       for(int x=0 ; x<workFieldDim.x ; ++x){
00671         //             concentrationArray[x][0]=concentrationArray[x][1];
00672         //             concentrationArray[x][workFieldDim.y-1]=concentrationArray[x][workFieldDim.y-2];
00673         //       }
00674         // 
00675         // 
00676         //    }
00677 
00678 
00679 }
00680 
00682 
00683 void FastDiffusionSolver2DFE::diffuseSingleField(unsigned int idx){
00684 
00685         // cerr<<" DIFFUSE"<<endl;
00686 
00687         // clock_t start, end;
00688         // double time;
00689         // start = clock();
00691 
00695 
00700 
00701 
00702         // Point3D pt, n;
00703         // unsigned int token = 0;
00704         // double distance;
00705         CellG * currentCellPtr=0,*nCell=0;
00706 
00707         // short currentCellType=0;
00708         float concentrationSum=0.0;
00709         float updatedConcentration=0.0;
00710 
00711 
00712         float currentConcentration=0.0;
00713         // short neighborCounter=0;
00714 
00715         DiffusionData & diffData = diffSecrFieldTuppleVec[idx].diffData;
00716         float diffConst=diffData.diffConst;
00717         float decayConst=diffData.decayConst;
00718         float deltaT=diffData.deltaT;
00719         float deltaX=diffData.deltaX;
00720         float dt_dx2=deltaT/(deltaX*deltaX);
00721 
00722 
00723 
00724         // std::set<unsigned char>::iterator sitr;
00725         // std::set<unsigned char>::iterator end_sitr=diffDataVec[idx].avoidTypeIdSet.end();
00726         // std::set<unsigned char>::iterator end_sitr_decay=diffDataVec[idx].avoidDecayInIdSet.end();
00727 
00728         // Automaton *automaton=potts->getAutomaton();
00729 
00730         ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
00731         ConcentrationField_t * scratchFieldPtr;
00732 
00733         //    scratchFieldPtr=concentrationFieldVector[diffDataVec.size()];
00734         if(!haveCouplingTerms)
00735                 scratchFieldPtr=concentrationFieldVector[diffSecrFieldTuppleVec.size()];
00736         else
00737                 scratchFieldPtr=concentrationFieldVector[diffSecrFieldTuppleVec.size()+idx];
00738 
00739         Array2D_t & concentrationArray = concentrationFieldPtr->getContainer();
00740         Array2D_t & scratchArray = scratchFieldPtr->getContainer();
00741         boundaryConditionInit(concentrationFieldPtr);//initializing boundary conditions
00742 
00743         //    cerr<<" concentrationArray[1][1] = "<<concentrationArray[1][1]<<endl;
00744         //    cerr<<"concentrationFieldPtr->get(Point3D(0,0,0))="<<concentrationFieldPtr->get(Point3D(0,0,0))<<endl;
00745 
00746         std::set<unsigned char>::iterator end_sitr_decay=diffData.avoidDecayInIdSet.end();   
00747         bool avoidDecayInMedium=false;
00748         //the assumption is that medium has type ID 0
00749 
00750         if(diffData.avoidDecayInIdSet.find(automaton->getTypeId("Medium")) != end_sitr_decay){
00751                 avoidDecayInMedium=true;
00752         }
00753 
00754 
00755 
00756         unsigned int offset=workFieldDim.x*workFieldDim.y;
00757         unsigned int counter=offset;
00758         float *temp;
00759 
00760         unsigned x_min=1,x_max=workFieldDim.x-1;
00761         unsigned y_min=1,y_max=workFieldDim.y-1;
00762 
00763         if(diffData.useBoxWatcher){
00764                 Point3D minCoordinates=*(boxWatcherSteppable->getMinCoordinatesPtr());
00765                 Point3D maxCoordinates=*(boxWatcherSteppable->getMaxCoordinatesPtr());
00766 
00767                 x_min=minCoordinates.x+1;
00768                 x_max=maxCoordinates.x+1;
00769                 y_min=minCoordinates.y+1;
00770                 y_max=maxCoordinates.y+1;
00771         }
00772 
00773         //       for (int x = 1; x < workFieldDim.x-1; ++x)
00774         //          for (int y = 1; y < workFieldDim.y-1; ++y){
00775         for (int x = x_min; x < x_max; ++x)
00776                 for (int y = y_min; y < y_max; ++y){
00777 
00778                         Point3D pt(x-1,y-1,0);
00779                         currentCellPtr=cellFieldG->get(pt);
00780                         currentConcentration = concentrationArray[x][y];
00781                         temp=&scratchArray[x][y];
00782                         //          temp  = &scratchVec[x][y][1];
00783                         //          temp=&scratch[counter];
00784                         //          scratchArray[x][y][1]=0;
00785 
00786                         updatedConcentration=0.0;
00787                         concentrationSum=0.0;
00788 
00789                         concentrationSum+=concentrationArray[x][y+1];
00790                         concentrationSum+=concentrationArray[x+1][y];
00791                         concentrationSum+=concentrationArray[x][y-1];
00792                         concentrationSum+=concentrationArray[x-1][y];
00793 
00794                         updatedConcentration =  dt_dx2*diffConst*(concentrationSum - 4*currentConcentration)+currentConcentration;
00795 
00796 
00797 
00798                         //processing decay depandent on type of the current cell
00799 
00800                         //          updatedConcentration-=deltaT*(decayConst*currentConcentration);//decay in this type is allowed
00801                         if(currentCellPtr){
00802                                 if(diffData.avoidDecayInIdSet.find(currentCellPtr->type)!=end_sitr_decay){
00803                                         ;//decay in this type is forbidden
00804                                 }else{
00805                                         updatedConcentration-=deltaT*(decayConst*currentConcentration);//decay in this type is allowed
00806                                 }
00807                         }else{
00808                                 if(avoidDecayInMedium){
00809                                         ;//decay in Medium is forbidden
00810                                 }else{
00811                                         updatedConcentration-=deltaT*(decayConst*currentConcentration); //decay in Medium is allowed
00812                                 }
00813                         }
00814                         //          if(!currentCellPtr){
00815                         //             updatedConcentration-=deltaT*(decayConst*currentConcentration);//decay in this type is allowed  
00816                         //          }
00817                         //          if (x<=3&&y<=3){
00818                         //             cerr<<"concentrationSum="<<concentrationSum<<" updatedConcentration="<<updatedConcentration<<endl;
00819                         //          }
00820 
00821                         *temp=updatedConcentration;
00822                         //               scratch[counter]=updatedConcentration;
00823                         //             ++counter;
00824                         //          scratchArray[x][y][1]=updatedConcentration;//updating scratch
00825 
00826                 }
00827 
00828 
00829 }
00830 
00832 void FastDiffusionSolver2DFE::diffuse() {
00833 
00834         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){
00835                 diffuseSingleField(i);
00836                 if(!haveCouplingTerms){ //without coupling terms field swap takes place immediately aftera given field has been diffused
00837                         ConcentrationField_t * concentrationField=concentrationFieldVector[i];
00838                         ConcentrationField_t * scratchField=concentrationFieldVector[diffSecrFieldTuppleVec.size()];
00839                         //copy updated values from scratch to concentration field
00840                         scrarch2Concentration(scratchField,concentrationField);
00841                 }
00842 
00843         }
00844 
00845         if(haveCouplingTerms){  //with coupling terms we swap scratch and concentration field after all fields have been diffused
00846                 for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){
00847                         ConcentrationField_t * concentrationField=concentrationFieldVector[i];
00848                         ConcentrationField_t * scratchField=concentrationFieldVector[diffSecrFieldTuppleVec.size()+i];
00849                         scrarch2Concentration(scratchField,concentrationField);
00850 
00851                 }
00852 
00853         }
00854 
00855 
00856 }
00858 void FastDiffusionSolver2DFE::scrarch2Concentration( ConcentrationField_t *scratchField, ConcentrationField_t *concentrationField){
00859         scratchField->switchContainersQuick(*(concentrationField));
00860 
00861 }
00862 
00864 void FastDiffusionSolver2DFE::outputField( std::ostream & _out, ConcentrationField_t *_concentrationField){
00865         Point3D pt;
00866         float tempValue;
00867 
00868 
00869         for (pt.x = 0; pt.x < fieldDim.x; pt.x++)
00870                 for (pt.y = 0; pt.y < fieldDim.y; pt.y++)
00871                         for (pt.z = 0; pt.z < fieldDim.z; pt.z++){
00872                                 tempValue=_concentrationField->get(pt);
00873                                 _out<<pt.x<<" "<<pt.y<<" "<<pt.z<<" "<<tempValue<<endl;
00874                         }
00875 }
00877 void FastDiffusionSolver2DFE::readConcentrationField(std::string fileName,ConcentrationField_t *concentrationField){
00878 
00879         ifstream in(fileName.c_str());
00880 
00881         ASSERT_OR_THROW(string("Could not open chemical concentration file '") +
00882                 fileName + "'!", in.is_open());
00883 
00884         Point3D pt;
00885         float c;
00886         //Zero entire field
00887         for (pt.z = 0; pt.z < fieldDim.z; pt.z++)
00888                 for (pt.y = 0; pt.y < fieldDim.y; pt.y++)
00889                         for (pt.x = 0; pt.x < fieldDim.x; pt.x++){
00890                                 concentrationField->set(pt,0);
00891                         }  
00892 
00893                         while(!in.eof()){
00894                                 in>>pt.x>>pt.y>>pt.z>>c;
00895                                 if(!in.fail())
00896                                         concentrationField->set(pt,c);
00897                         }
00898 
00899 }
00900 
00901 
00902 void FastDiffusionSolver2DFE::update(CC3DXMLElement *_xmlData, bool _fullInitFlag){
00903 
00904         //notice, only basic steering is enabled for PDE solvers - changing diffusion constants, do -not-diffuse to types etc...
00905         // Coupling coefficients cannot be changed and also there is no way to allocate extra fields while simulation is running
00906         diffSecrFieldTuppleVec.clear();
00907 
00908         CC3DXMLElementList diffFieldXMLVec=_xmlData->getElements("DiffusionField");
00909         for(unsigned int i = 0 ; i < diffFieldXMLVec.size() ; ++i ){
00910                 diffSecrFieldTuppleVec.push_back(DiffusionSecretionFastFieldTupple());
00911                 DiffusionData & diffData=diffSecrFieldTuppleVec[diffSecrFieldTuppleVec.size()-1].diffData;
00912                 SecretionData & secrData=diffSecrFieldTuppleVec[diffSecrFieldTuppleVec.size()-1].secrData;
00913 
00914                 if(diffFieldXMLVec[i]->findElement("DiffusionData"))
00915                         diffData.update(diffFieldXMLVec[i]->getFirstElement("DiffusionData"));
00916 
00917                 if(diffFieldXMLVec[i]->findElement("SecretionData"))
00918                         secrData.update(diffFieldXMLVec[i]->getFirstElement("SecretionData"));
00919 
00920                 if(diffFieldXMLVec[i]->findElement("ReadFromFile"))
00921                         readFromFileFlag=true;
00922 
00923 
00924         }
00925         if(_xmlData->findElement("Serialize")){
00926 
00927                 serializeFlag=true;
00928                 if(_xmlData->getFirstElement("Serialize")->findAttribute("Frequency")){
00929                         serializeFrequency=_xmlData->getFirstElement("Serialize")->getAttributeAsUInt("Frequency");
00930                 }
00931                 cerr<<"serialize Flag="<<serializeFlag<<endl;
00932 
00933         }
00934 
00935         if(_xmlData->findElement("ReadFromFile")){
00936                 readFromFileFlag=true;
00937                 cerr<<"readFromFileFlag="<<readFromFileFlag<<endl;
00938         }
00939 
00940 
00941 
00942 
00943 
00944         for(int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00945                 diffSecrFieldTuppleVec[i].diffData.setAutomaton(automaton);
00946                 diffSecrFieldTuppleVec[i].secrData.setAutomaton(automaton);
00947                 diffSecrFieldTuppleVec[i].diffData.initialize(automaton);
00948                 diffSecrFieldTuppleVec[i].secrData.initialize(automaton);
00949         }
00950 
00952 
00953         for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
00954 
00955                 diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec.assign(diffSecrFieldTuppleVec[i].secrData.secrTypesNameSet.size(),0);
00956                 unsigned int j=0;
00957                 for(set<string>::iterator sitr=diffSecrFieldTuppleVec[i].secrData.secrTypesNameSet.begin() ; sitr != diffSecrFieldTuppleVec[i].secrData.secrTypesNameSet.end()  ; ++sitr){
00958 
00959                         if((*sitr)=="Secretion"){
00960                                 diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec[j]=&FastDiffusionSolver2DFE::secreteSingleField;
00961                                 ++j;
00962                         }
00963                         else if((*sitr)=="SecretionOnContact"){
00964                                 diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec[j]=&FastDiffusionSolver2DFE::secreteOnContactSingleField;
00965                                 ++j;
00966                         }
00967                         else if((*sitr)==