From a56de07145a83fda6241e1c6ff438fa3e1c8837d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Asia=20J=C4=99drzejewska-Szmek?= <asia@in.waw.pl> Date: Fri, 20 Jan 2017 15:49:30 -0500 Subject: [PATCH] Fix bugs: DifShell now compiles and seems to give correct results --- moose-core/biophysics/DifBuffer.cpp | 164 ++++++++++++- moose-core/biophysics/DifBuffer.h | 89 +++---- moose-core/biophysics/DifBufferBase.cpp | 178 ++++++++------ moose-core/biophysics/DifBufferBase.h | 46 +++- moose-core/biophysics/DifShell.cpp | 133 +++++++++- moose-core/biophysics/DifShell.h | 30 ++- moose-core/biophysics/DifShellBase.cpp | 311 ++++++++++-------------- moose-core/biophysics/DifShellBase.h | 37 ++- 8 files changed, 649 insertions(+), 339 deletions(-) diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp index 1d6ccf9d..0d695f80 100644 --- a/moose-core/biophysics/DifBuffer.cpp +++ b/moose-core/biophysics/DifBuffer.cpp @@ -46,6 +46,7 @@ // Code: #include "header.h" +#include "DifBufferBase.h" #include "DifBuffer.h" #include "ElementValueFinfo.h" #include "../utility/numutil.h" @@ -65,7 +66,7 @@ const Cinfo * DifBuffer::initCinfo() static Dinfo<DifBuffer> dinfo; static Cinfo difBufferCinfo( "DifBuffer", - Neutral::initCinfo(), + DifBufferBase::initCinfo(), 0, 0, &dinfo, @@ -83,13 +84,23 @@ static const Cinfo * difBufferCinfo = DifBuffer::initCinfo(); //////////////////////////////////////////////////////////////////////////////// DifBuffer::DifBuffer() : + activation_(0), + Af_(0), + Bf_(0), + bFree_(0), + bBound_(0), + bTot_(0), + kf_(0), + kb_(0), + D_(0), shapeMode_(0), - diameter_(0), length_(0), + diameter_(0), thickness_(0), + volume_(0), outerArea_(0), - innerArea_(0), - volume_(0) + innerArea_(0) + {} //////////////////////////////////////////////////////////////////////////////// @@ -114,12 +125,40 @@ double DifBuffer::vGetBFree(const Eref& e) const { return bFree_; } +void DifBuffer::vSetBFree(const Eref& e,double value) +{ + if ( value < 0.0 ) { + cerr << "Error: DifBuffer: Free Buffer cannot be negative!\n"; + return; + } + if (value > bTot_){ + cerr << "Error: DifBuffer: Free Buffer cannot exceed total buffer!\n"; + return; + } + bFree_ = value; + bBound_ = bTot_ - bFree_; +} double DifBuffer::vGetBBound(const Eref& e) const { return bBound_; } +void DifBuffer::vSetBBound(const Eref& e,double value) +{ + if ( value < 0.0 ) { + cerr << "Error: DifBuffer: Bound Buffer cannot be negative!\n"; + return; + } + if (value > bTot_){ + cerr << "Error: DifBuffer: Bound Buffer cannot exceed total buffer!\n"; + return; + } + bBound_ = value; + bFree_ = bTot_ - bBound_; +} + + double DifBuffer::vGetBTot(const Eref& e) const { return bTot_; @@ -181,12 +220,129 @@ void DifBuffer::vSetD(const Eref& e,double value) } +void DifBuffer::vSetShapeMode(const Eref& e, unsigned int shapeMode ) +{ + if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) { + cerr << "Error: DifBuffer: I only understand shapeModes 0, 1 and 3.\n"; + return; + } + shapeMode_ = shapeMode; +} + +unsigned int DifBuffer::vGetShapeMode(const Eref& e) const +{ + return shapeMode_; +} + +void DifBuffer::vSetLength(const Eref& e, double length ) +{ + if ( length < 0.0 ) { + cerr << "Error: DifBuffer: length cannot be negative!\n"; + return; + } + + length_ = length; +} + +double DifBuffer::vGetLength(const Eref& e ) const +{ + return length_; +} + +void DifBuffer::vSetDiameter(const Eref& e, double diameter ) +{ + if ( diameter < 0.0 ) { + cerr << "Error: DifBuffer: diameter cannot be negative!\n"; + return; + } + + diameter_ = diameter; +} + +double DifBuffer::vGetDiameter(const Eref& e ) const +{ + return diameter_; +} + +void DifBuffer::vSetThickness( const Eref& e, double thickness ) +{ + if ( thickness < 0.0 ) { + cerr << "Error: DifBuffer: thickness cannot be negative!\n"; + return; + } + + thickness_ = thickness; +} + +double DifBuffer::vGetThickness(const Eref& e) const +{ + return thickness_; +} + +void DifBuffer::vSetVolume(const Eref& e, double volume ) +{ + if ( shapeMode_ != 3 ) + cerr << "Warning: DifBuffer: Trying to set volume, when shapeMode is not USER-DEFINED\n"; + + if ( volume < 0.0 ) { + cerr << "Error: DifBuffer: volume cannot be negative!\n"; + return; + } + + volume_ = volume; +} + +double DifBuffer::vGetVolume(const Eref& e ) const +{ + return volume_; +} + +void DifBuffer::vSetOuterArea(const Eref& e, double outerArea ) +{ + if (shapeMode_ != 3 ) + cerr << "Warning: DifBuffer: Trying to set outerArea, when shapeMode is not USER-DEFINED\n"; + + if ( outerArea < 0.0 ) { + cerr << "Error: DifBuffer: outerArea cannot be negative!\n"; + return; + } + + outerArea_ = outerArea; +} + +double DifBuffer::vGetOuterArea(const Eref& e ) const +{ + return outerArea_; +} + +void DifBuffer::vSetInnerArea(const Eref& e, double innerArea ) +{ + if ( shapeMode_ != 3 ) + cerr << "Warning: DifBuffer: Trying to set innerArea, when shapeMode is not USER-DEFINED\n"; + + if ( innerArea < 0.0 ) { + cerr << "Error: DifBuffer: innerArea cannot be negative!\n"; + return; + } + + innerArea_ = innerArea; +} + +double DifBuffer::vGetInnerArea(const Eref& e) const +{ + return innerArea_; +} + + + + void DifBuffer::vBuffer(const Eref& e, double C ) { activation_ = C; } + double DifBuffer::integrate( double state, double dt, double A, double B ) { if ( B > EPSILON ) { diff --git a/moose-core/biophysics/DifBuffer.h b/moose-core/biophysics/DifBuffer.h index 57307491..fe3f7772 100644 --- a/moose-core/biophysics/DifBuffer.h +++ b/moose-core/biophysics/DifBuffer.h @@ -1,54 +1,16 @@ -/* DifBuffer.h --- - * - * Filename: DifBuffer.h - * Description: - * Author: Subhasis Ray - * Maintainer: - * Created: Mon Feb 16 11:42:50 2015 (-0500) - * Version: - * Package-Requires: () - * Last-Updated: Mon Feb 23 10:49:54 2015 (-0500) - * By: Subhasis Ray - * Update #: 25 - * URL: - * Doc URL: - * Keywords: - * Compatibility: - * - */ - -/* Commentary: - * - * - * - */ - -/* Change Log: - * - * - */ - -/* This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or (at - * your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>. - */ - -/* Code: */ - +/********************************************************************** + ** This program is part of 'MOOSE', the + ** Multiscale Object Oriented Simulation Environment. + ** copyright (C) 2003-2008 + ** Upinder S. Bhalla, Niraj Dudani and NCBS + ** It is made available under the terms of the + ** GNU Lesser General Public License version 2.1 + ** See the file COPYING.LIB for the full notice. + *****/ #ifndef _DifBuffer_h #define _DifBuffer_h -class DifBuffer -{ +class DifBuffer: public DifBufferBase{ public: DifBuffer(); @@ -63,7 +25,10 @@ class DifBuffer void vSetActivation(const Eref& e,double value); double vGetBFree(const Eref& e) const; + void vSetBFree(const Eref& e,double value); + double vGetBBound(const Eref& e) const; + void vSetBBound(const Eref& e,double value); double vGetBTot(const Eref& e) const; // total buffer concentration in mM (free + bound) void vSetBTot(const Eref& e,double value); @@ -76,10 +41,27 @@ class DifBuffer double vGetD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec void vSetD(const Eref& e,double value); + void vSetShapeMode(const Eref& e, unsigned int shapeMode ); + unsigned int vGetShapeMode(const Eref& e) const; + + void vSetLength(const Eref& e, double length ); + double vGetLength(const Eref& e) const; + + void vSetDiameter(const Eref& e, double diameter ); + double vGetDiameter(const Eref& e) const; + + void vSetThickness(const Eref& e, double thickness ); + double vGetThickness(const Eref& e) const; - //void concentration(); + void vSetVolume(const Eref& e, double volume ); + double vGetVolume(const Eref& e) const; + + void vSetOuterArea(const Eref& e, double outerArea ); + double vGetOuterArea(const Eref& e) const; - + void vSetInnerArea(const Eref& e, double innerArea ); + double vGetInnerArea(const Eref& e) const; + static const Cinfo * initCinfo(); private: @@ -98,6 +80,13 @@ class DifBuffer double kf_; //forward rate constant double kb_; //backward rate constant double D_; //diffusion constant + unsigned int shapeMode_; + double length_; + double diameter_; + double thickness_; + double volume_; + double outerArea_; + double innerArea_; static const double EPSILON; }; diff --git a/moose-core/biophysics/DifBufferBase.cpp b/moose-core/biophysics/DifBufferBase.cpp index 3b894c35..72725ac4 100644 --- a/moose-core/biophysics/DifBufferBase.cpp +++ b/moose-core/biophysics/DifBufferBase.cpp @@ -137,7 +137,7 @@ const Cinfo * DifBufferBase::initCinfo() "Diameter of shell", &DifBufferBase::setDiameter, &DifBufferBase::getDiameter); - static ElementValueFinfo<DifBufferBase, int> shapeMode("shapeMode", + static ElementValueFinfo<DifBufferBase, unsigned int> shapeMode("shapeMode", "shape of the shell: SHELL=0, SLICE=SLAB=1, USERDEF=3", &DifBufferBase::setShapeMode, &DifBufferBase::getShapeMode); @@ -191,9 +191,9 @@ const Cinfo * DifBufferBase::initCinfo() &innerDif, &outerDif, // - //reactionOut(), - //innerDifSourceOut(), - //outerDifSourceOut(), + reactionOut(), + innerDifSourceOut(), + outerDifSourceOut(), ////////////////////////////////////////////////////////////////// // DestFinfo definitions ////////////////////////////////////////////////////////////////// @@ -206,7 +206,7 @@ const Cinfo * DifBufferBase::initCinfo() "Description", "Models diffusible buffer where total concentration is constant. It is" " coupled with a DifShell.", }; - static Dinfo<DifBufferBase> dinfo; + static ZeroSizeDinfo<int> dinfo; static Cinfo difBufferCinfo( "DifBufferBase", Neutral::initCinfo(), @@ -226,15 +226,8 @@ static const Cinfo * difBufferCinfo = DifBufferBase::initCinfo(); // Class functions //////////////////////////////////////////////////////////////////////////////// -DifBufferBase::DifBufferBase() : - shapeMode_(0), - diameter_(0), - length_(0), - thickness_(0), - outerArea_(0), - innerArea_(0), - volume_(0) -{} +DifBufferBase::DifBufferBase() +{ ; } double DifBufferBase::getActivation(const Eref& e) const @@ -253,10 +246,19 @@ double DifBufferBase::getBFree(const Eref& e) const return vGetBFree(e); } +void DifBufferBase::setBFree(const Eref& e,double value) +{ + vSetBFree(e,value); +} + double DifBufferBase::getBBound(const Eref& e) const { return vGetBBound(e); } +void DifBufferBase::setBBound(const Eref& e,double value) +{ + vSetBBound(e,value); +} double DifBufferBase::getBTot(const Eref& e) const { @@ -299,117 +301,77 @@ void DifBufferBase::setD(const Eref& e,double value) vSetD(e,value); } -int DifBufferBase::getShapeMode(const Eref& e) const +void DifBufferBase::setShapeMode(const Eref& e, unsigned int shapeMode ) { - return shapeMode_; + vSetShapeMode(e,shapeMode); } -void DifBufferBase::setShapeMode(const Eref& e,int value) +unsigned int DifBufferBase::getShapeMode(const Eref& e) const { - if ( value != 0 && value !=1 && value != 3 ) { - cerr << "Error: DifBuffer: Shape mode can only be 0, 1 or 3"; - return; - } - shapeMode_ = value; + return vGetShapeMode(e); } -double DifBufferBase::getLength(const Eref& e) const +void DifBufferBase::setLength(const Eref& e, double length ) { - return length_; + vSetLength(e,length); } -void DifBufferBase::setLength(const Eref& e,double value) +double DifBufferBase::getLength(const Eref& e ) const { - if ( value < 0.0) { - cerr << "Error: DifBuffer: Length cannot be negative!\n"; - return; - } - length_ = value; + return vGetLength(e); } -double DifBufferBase::getDiameter(const Eref& e) const +void DifBufferBase::setDiameter(const Eref& e, double diameter ) { - return diameter_; + vSetDiameter(e,diameter); } -void DifBufferBase::setDiameter(const Eref& e,double value) +double DifBufferBase::getDiameter(const Eref& e ) const { - if ( value < 0.0) { - cerr << "Error: DifBuffer: Diameter cannot be negative!\n"; - return; - } - diameter_ = value; + return vGetDiameter(e); } -double DifBufferBase::getThickness(const Eref& e) const +void DifBufferBase::setThickness( const Eref& e, double thickness ) { - return thickness_; + vSetThickness(e,thickness); } -void DifBufferBase::setThickness(const Eref& e,double value) +double DifBufferBase::getThickness(const Eref& e) const { - if ( value < 0.0) { - cerr << "Error: DifBuffer: Thickness cannot be negative!\n"; - return; - } - thickness_ = value; + return vGetThickness(e); } void DifBufferBase::setVolume(const Eref& e, double volume ) { - if ( shapeMode_ != 3 ) - cerr << "Warning: DifBuffer: Trying to set volume, when shapeMode is not USER-DEFINED\n"; - - if ( volume < 0.0 ) { - cerr << "Error: DifBuffer: volume cannot be negative!\n"; - return; - } - - volume_ = volume; + vSetVolume(e,volume); } double DifBufferBase::getVolume(const Eref& e ) const { - return volume_; + return vGetVolume(e); } void DifBufferBase::setOuterArea(const Eref& e, double outerArea ) { - if (shapeMode_ != 3 ) - cerr << "Warning: DifBuffer: Trying to set outerArea, when shapeMode is not USER-DEFINED\n"; - - if ( outerArea < 0.0 ) { - cerr << "Error: DifBuffer: outerArea cannot be negative!\n"; - return; - } - - outerArea_ = outerArea; + vSetOuterArea(e,outerArea); } double DifBufferBase::getOuterArea(const Eref& e ) const { - return outerArea_; + return vGetOuterArea(e); } void DifBufferBase::setInnerArea(const Eref& e, double innerArea ) { - if ( shapeMode_ != 3 ) - cerr << "Warning: DifBuffer: Trying to set innerArea, when shapeMode is not USER-DEFINED\n"; - - if ( innerArea < 0.0 ) { - cerr << "Error: DifBuffer: innerArea cannot be negative!\n"; - return; - } - - innerArea_ = innerArea; + vSetInnerArea(e,innerArea); } double DifBufferBase::getInnerArea(const Eref& e) const { - return innerArea_; + return vGetInnerArea(e); } -} + void DifBufferBase::buffer(const Eref& e,double C) { @@ -433,3 +395,65 @@ void DifBufferBase:: fluxFromIn(const Eref& e,double innerC, double innerThickne { vFluxFromIn(e,innerC,innerThickness); } +void DifBufferBase::vSetSolver( const Eref& e, Id hsolve ) +{;} + +void DifBufferBase::zombify( Element* orig, const Cinfo* zClass, + Id hsolve ) +{ + if ( orig->cinfo() == zClass ) + return; + unsigned int start = orig->localDataStart(); + unsigned int num = orig->numLocalData(); + if ( num == 0 ) + return; + unsigned int len = 14; + vector< double > data( num * len ); + + unsigned int j = 0; + + for ( unsigned int i = 0; i < num; ++i ) { + Eref er( orig, i + start ); + const DifBufferBase* ds = + reinterpret_cast< const DifBufferBase* >( er.data() ); + data[j + 0] = ds->getActivation( er ); + data[j + 1] = ds->getBFree( er ); + data[j + 2] = ds->getBBound( er ); + data[j + 3] = ds->getBTot( er ); + data[j + 4] = ds->getKf( er ); + data[j + 5] = ds->getKb( er ); + data[j + 6] = ds->getD( er ); + data[j + 7] = ds->getShapeMode( er ); + data[j + 8] = ds->getLength( er ); + data[j + 9] = ds->getDiameter( er ); + data[j + 10] = ds->getThickness( er ); + data[j + 11] = ds->getVolume( er ); + data[j + 12] = ds->getOuterArea( er ); + data[j + 13] = ds->getInnerArea( er ); + j += len; + } + orig->zombieSwap( zClass ); + j = 0; + for ( unsigned int i = 0; i < num; ++i ) { + Eref er( orig, i + start ); + DifBufferBase* ds = + reinterpret_cast< DifBufferBase* >( er.data() ); + ds->vSetSolver(er,hsolve); + ds->setActivation(er, data[j+0]); + ds->setBFree(er, data[j + 1]); + ds->setBBound(er, data[j + 2]); + ds->setBTot(er, data[j + 3]); + ds->setKf(er, data[j + 4]); + ds->setKb(er, data[j + 5]); + ds->setD(er, data[j + 6]); + ds->setShapeMode(er, data[j + 7]); + ds->setLength(er, data[j + 8]); + ds->setDiameter(er, data[j + 9]); + ds->setThickness(er, data[j + 10]); + ds->setVolume(er, data[j + 11]); + ds->setOuterArea(er, data[j + 12]); + ds->setInnerArea(er, data[j + 13]); + j += len; //?? + } + +} diff --git a/moose-core/biophysics/DifBufferBase.h b/moose-core/biophysics/DifBufferBase.h index e169ec29..e1bc23b6 100644 --- a/moose-core/biophysics/DifBufferBase.h +++ b/moose-core/biophysics/DifBufferBase.h @@ -14,6 +14,7 @@ class DifBufferBase { public: + DifBufferBase(); void buffer(const Eref& e,double C); void reinit( const Eref & e, ProcPtr p ); @@ -31,8 +32,11 @@ public: void setActivation(const Eref& e,double value); double getBFree(const Eref& e) const; + void setBFree(const Eref& e,double value); + double getBBound(const Eref& e) const; - + void setBBound(const Eref& e,double value); + double getBTot(const Eref& e) const; // total buffer concentration in mM (free + bound) void setBTot(const Eref& e,double value); @@ -46,8 +50,8 @@ public: void setD(const Eref& e,double value); - int getShapeMode(const Eref& e) const; - void setShapeMode(const Eref& e,int value); // variables SHELL=0, SLICE=SLAB=1, USERDEF=3. + unsigned int getShapeMode(const Eref& e) const; + void setShapeMode(const Eref& e,unsigned int value); // variables SHELL=0, SLICE=SLAB=1, USERDEF=3. double getLength(const Eref& e) const; // shell length void setLength(const Eref& e,double value); @@ -72,7 +76,11 @@ public: virtual void vSetActivation(const Eref& e,double value) = 0; virtual double vGetBFree(const Eref& e) const = 0; + virtual void vSetBFree(const Eref& e,double value) = 0; + virtual double vGetBBound(const Eref& e) const = 0; + virtual void vSetBBound(const Eref& e,double value) = 0; + virtual double vGetBTot(const Eref& e) const = 0; // total buffer concentration in mM (free + bound) virtual void vSetBTot(const Eref& e,double value) = 0; @@ -86,19 +94,37 @@ public: virtual double vGetD(const Eref& e) const = 0; // diffusion constant of buffer molecules, m^2/sec virtual void vSetD(const Eref& e,double value) = 0; + virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0; + virtual unsigned int vGetShapeMode(const Eref& e) const = 0; + + virtual void vSetLength(const Eref& e, double length ) = 0; + virtual double vGetLength(const Eref& e) const = 0; + + virtual void vSetDiameter(const Eref& e, double diameter ) = 0; + virtual double vGetDiameter(const Eref& e) const = 0; + + virtual void vSetThickness(const Eref& e, double thickness ) = 0; + virtual double vGetThickness(const Eref& e) const = 0; + + virtual void vSetVolume(const Eref& e, double volume ) = 0; + virtual double vGetVolume(const Eref& e) const = 0; + + virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0; + virtual double vGetOuterArea(const Eref& e) const = 0; + + virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0; + virtual double vGetInnerArea(const Eref& e) const = 0; + static SrcFinfo4< double, double, double, double >* reactionOut(); static SrcFinfo2< double, double >* innerDifSourceOut(); static SrcFinfo2< double, double >* outerDifSourceOut(); static const Cinfo * initCinfo(); + virtual void vSetSolver( const Eref& e, Id hsolve ); + static void zombify( Element* orig, const Cinfo* zClass, + Id hsolve ); private: - unsigned int shapeMode_; - double length_; - double diameter_; - double thickness_; - double volume_; - double outerArea_; - double innerArea_; + }; diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp index b029f3a1..a9f9ee77 100644 --- a/moose-core/biophysics/DifShell.cpp +++ b/moose-core/biophysics/DifShell.cpp @@ -9,8 +9,9 @@ **********************************************************************/ #include "header.h" -#include "DifShell.h" #include "DifShellBase.h" +#include "DifShell.h" + const double DifShell::EPSILON = 1.0e-10; const double DifShell::F = 96485.3415; /* C / mol like in genesis */ @@ -32,7 +33,7 @@ const Cinfo* DifShell::initCinfo() static Dinfo< DifShell > dinfo; static Cinfo difShellCinfo( "DifShell", - Neutral::initCinfo(), + DifShellBase::initCinfo(), 0, 0, &dinfo, @@ -59,7 +60,14 @@ DifShell::DifShell() : Ceq_( 0.0 ), D_( 0.0 ), valence_( 0.0 ), - leak_( 0.0 ) + leak_( 0.0 ), + shapeMode_(0), + length_(0), + diameter_(0), + thickness_(0), + volume_(0), + outerArea_(0), + innerArea_(0) { ; } //////////////////////////////////////////////////////////////////////////////// @@ -137,6 +145,122 @@ double DifShell::vGetLeak(const Eref& e ) const return leak_; } + +void DifShell::vSetShapeMode(const Eref& e, unsigned int shapeMode ) +{ + if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) { + cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n"; + return; + } + shapeMode_ = shapeMode; +} + +unsigned int DifShell::vGetShapeMode(const Eref& e) const +{ + return shapeMode_; +} + +void DifShell::vSetLength(const Eref& e, double length ) +{ + if ( length < 0.0 ) { + cerr << "Error: DifShell: length cannot be negative!\n"; + return; + } + + length_ = length; +} + +double DifShell::vGetLength(const Eref& e ) const +{ + return length_; +} + +void DifShell::vSetDiameter(const Eref& e, double diameter ) +{ + if ( diameter < 0.0 ) { + cerr << "Error: DifShell: diameter cannot be negative!\n"; + return; + } + + diameter_ = diameter; +} + +double DifShell::vGetDiameter(const Eref& e ) const +{ + return diameter_; +} + +void DifShell::vSetThickness( const Eref& e, double thickness ) +{ + if ( thickness < 0.0 ) { + cerr << "Error: DifShell: thickness cannot be negative!\n"; + return; + } + + thickness_ = thickness; +} + +double DifShell::vGetThickness(const Eref& e) const +{ + return thickness_; +} + +void DifShell::vSetVolume(const Eref& e, double volume ) +{ + if ( shapeMode_ != 3 ) + cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n"; + + if ( volume < 0.0 ) { + cerr << "Error: DifShell: volume cannot be negative!\n"; + return; + } + + volume_ = volume; +} + +double DifShell::vGetVolume(const Eref& e ) const +{ + return volume_; +} + +void DifShell::vSetOuterArea(const Eref& e, double outerArea ) +{ + if (shapeMode_ != 3 ) + cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n"; + + if ( outerArea < 0.0 ) { + cerr << "Error: DifShell: outerArea cannot be negative!\n"; + return; + } + + outerArea_ = outerArea; +} + +double DifShell::vGetOuterArea(const Eref& e ) const +{ + return outerArea_; +} + +void DifShell::vSetInnerArea(const Eref& e, double innerArea ) +{ + if ( shapeMode_ != 3 ) + cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n"; + + if ( innerArea < 0.0 ) { + cerr << "Error: DifShell: innerArea cannot be negative!\n"; + return; + } + + innerArea_ = innerArea; +} + +double DifShell::vGetInnerArea(const Eref& e) const +{ + return innerArea_; +} + + + //////////////////////////////////////////////////////////////////////////////// // Local dest functions //////////////////////////////////////////////////////////////////////////////// @@ -236,6 +360,7 @@ void DifShell::vFluxFromOut(const Eref& e, double outerC, double outerThickness * We could pre-compute ( D / Volume ), but let us leave the optimizations * for the solver. */ + //cout << "FluxFromOut "<<outerC<<" "<<outerThickness; dCbyDt_ += diff * outerC; Cmultiplier_ += diff ; } @@ -245,7 +370,7 @@ void DifShell::vFluxFromIn(const Eref& e, double innerC, double innerThickness ) //influx from inner shell //double dx = ( innerThickness + thickness_ ) / 2.0; double diff = 2.*( D_ / volume_ ) * ( innerArea_ / (innerThickness + thickness_) ); - + //cout << "FluxFromIn "<<innerC<<" "<<innerThickness; dCbyDt_ += diff * innerC ; Cmultiplier_ += diff ; } diff --git a/moose-core/biophysics/DifShell.h b/moose-core/biophysics/DifShell.h index 3bc4d07c..6bef0803 100644 --- a/moose-core/biophysics/DifShell.h +++ b/moose-core/biophysics/DifShell.h @@ -11,8 +11,7 @@ #ifndef _DIFSHELL_H #define _DIFSHELL_H -class DifShell: public DifShellBase -{ +class DifShell: public DifShellBase{ public: DifShell(); @@ -54,6 +53,26 @@ class DifShell: public DifShellBase void vSetLeak(const Eref& e, double leak ); double vGetLeak(const Eref& e) const; + void vSetShapeMode(const Eref& e, unsigned int shapeMode ); + unsigned int vGetShapeMode(const Eref& e) const; + + void vSetLength(const Eref& e, double length ); + double vGetLength(const Eref& e) const; + + void vSetDiameter(const Eref& e, double diameter ); + double vGetDiameter(const Eref& e) const; + + void vSetThickness(const Eref& e, double thickness ); + double vGetThickness(const Eref& e) const; + + void vSetVolume(const Eref& e, double volume ); + double vGetVolume(const Eref& e) const; + + void vSetOuterArea(const Eref& e, double outerArea ); + double vGetOuterArea(const Eref& e) const; + + void vSetInnerArea(const Eref& e, double innerArea ); + double vGetInnerArea(const Eref& e) const; static const Cinfo * initCinfo(); @@ -68,6 +87,13 @@ class DifShell: public DifShellBase double D_; double valence_; double leak_; + unsigned int shapeMode_; + double length_; + double diameter_; + double thickness_; + double volume_; + double outerArea_; + double innerArea_; static const double EPSILON; /// Faraday's constant (Coulomb / Mole) diff --git a/moose-core/biophysics/DifShellBase.cpp b/moose-core/biophysics/DifShellBase.cpp index e084b7e4..2400e392 100644 --- a/moose-core/biophysics/DifShellBase.cpp +++ b/moose-core/biophysics/DifShellBase.cpp @@ -100,42 +100,42 @@ const Cinfo* DifShellBase::initCinfo() sizeof( outerDifShared ) / sizeof( Finfo* )); static ElementValueFinfo< DifShellBase, double> C( "C", - "Concentration C",// is computed by the DifShell", - &DifShellBase::setC, - &DifShellBase::getC); + "Concentration C",// is computed by the DifShell", + &DifShellBase::setC, + &DifShellBase::getC); static ElementValueFinfo< DifShellBase, double> Ceq( "Ceq", "", - &DifShellBase::setCeq, - &DifShellBase::getCeq); + &DifShellBase::setCeq, + &DifShellBase::getCeq); static ElementValueFinfo< DifShellBase, double> D( "D", "", - &DifShellBase::setD, - &DifShellBase::getD); + &DifShellBase::setD, + &DifShellBase::getD); static ElementValueFinfo< DifShellBase, double> valence( "valence", "", - &DifShellBase::setValence, - &DifShellBase::getValence); + &DifShellBase::setValence, + &DifShellBase::getValence); static ElementValueFinfo< DifShellBase, double> leak( "leak", "", - &DifShellBase::setLeak, - &DifShellBase::getLeak); + &DifShellBase::setLeak, + &DifShellBase::getLeak); static ElementValueFinfo< DifShellBase, unsigned int> shapeMode( "shapeMode", "", - &DifShellBase::setShapeMode, - &DifShellBase::getShapeMode); + &DifShellBase::setShapeMode, + &DifShellBase::getShapeMode); static ElementValueFinfo< DifShellBase, double> length( "length", "", - &DifShellBase::setLength, - &DifShellBase::getLength); + &DifShellBase::setLength, + &DifShellBase::getLength); static ElementValueFinfo< DifShellBase, double> diameter( "diameter", "", - &DifShellBase::setDiameter, - &DifShellBase::getDiameter ); + &DifShellBase::setDiameter, + &DifShellBase::getDiameter ); static ElementValueFinfo< DifShellBase, double> thickness( "thickness", "", - &DifShellBase::setThickness, - &DifShellBase::getThickness ); + &DifShellBase::setThickness, + &DifShellBase::getThickness ); static ElementValueFinfo< DifShellBase, double> volume( "volume", "", - &DifShellBase::setVolume, - &DifShellBase::getVolume ); + &DifShellBase::setVolume, + &DifShellBase::getVolume ); static ElementValueFinfo< DifShellBase, double> outerArea( "outerArea", "", - &DifShellBase::setOuterArea, - &DifShellBase::getOuterArea); + &DifShellBase::setOuterArea, + &DifShellBase::getOuterArea); static ElementValueFinfo< DifShellBase, double> innerArea( "innerArea", "", - &DifShellBase::setInnerArea, - &DifShellBase::getInnerArea ); + &DifShellBase::setInnerArea, + &DifShellBase::getInnerArea ); static DestFinfo influx( "influx", "", @@ -183,9 +183,9 @@ const Cinfo* DifShellBase::initCinfo() ////////////////////////////////////////////////////////////////// &proc, &buffer, - // concentrationOut(), - //innerDifSourceOut(), - //outerDifSourceOut(), + concentrationOut(), + innerDifSourceOut(), + outerDifSourceOut(), &innerDif, &outerDif, ////////////////////////////////////////////////////////////////// @@ -213,15 +213,16 @@ const Cinfo* DifShellBase::initCinfo() "changes in concentration due to pumping, buffering and channel currents, by " "talking to the appropriate objects.", }; - static Dinfo< DifShellBase> dinfo; + static ZeroSizeDinfo< int > dinfo; + //static Dinfo< DifShellBase> dinfo; static Cinfo difShellBaseCinfo( - "DifShellBase", - Neutral::initCinfo(), - difShellBaseFinfos, - sizeof( difShellBaseFinfos ) / sizeof( Finfo* ), - &dinfo, - doc, - sizeof( doc ) / sizeof( string )); + "DifShellBase", + Neutral::initCinfo(), + difShellBaseFinfos, + sizeof( difShellBaseFinfos ) / sizeof( Finfo* ), + &dinfo, + doc, + sizeof( doc ) / sizeof( string )); return &difShellBaseCinfo; } @@ -229,14 +230,7 @@ const Cinfo* DifShellBase::initCinfo() //Cinfo *object* corresponding to the class. static const Cinfo* difShellBaseCinfo = DifShellBase::initCinfo(); -DifShellBase::DifShellBase() : - shapeMode_( 0 ), - length_( 0.0 ), - diameter_( 0.0 ), - thickness_( 0.0 ), - volume_( 0.0 ), - outerArea_( 0.0 ), - innerArea_( 0.0 ) +DifShellBase::DifShellBase() { ; } void DifShellBase::setC( const Eref& e, double C) @@ -250,7 +244,7 @@ double DifShellBase::getC(const Eref& e) const void DifShellBase::setCeq( const Eref& e, double Ceq ) { - vSetCeq(e,C); + vSetCeq(e,Ceq); } double DifShellBase::getCeq(const Eref& e ) const @@ -260,7 +254,7 @@ double DifShellBase::getCeq(const Eref& e ) const void DifShellBase::setD(const Eref& e, double D ) { - vSetD(e,D) + vSetD(e,D); } double DifShellBase::getD(const Eref& e ) const @@ -290,115 +284,72 @@ double DifShellBase::getLeak(const Eref& e ) const void DifShellBase::setShapeMode(const Eref& e, unsigned int shapeMode ) { - if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) { - cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n"; - return; - } - shapeMode_ = shapeMode; + vSetShapeMode(e,shapeMode); } unsigned int DifShellBase::getShapeMode(const Eref& e) const { - return shapeMode_; + return vGetShapeMode(e); } void DifShellBase::setLength(const Eref& e, double length ) { - if ( length < 0.0 ) { - cerr << "Error: DifShell: length cannot be negative!\n"; - return; - } - - length_ = length; + vSetLength(e,length); } double DifShellBase::getLength(const Eref& e ) const { - return length_; + return vGetLength(e); } void DifShellBase::setDiameter(const Eref& e, double diameter ) { - if ( diameter < 0.0 ) { - cerr << "Error: DifShell: diameter cannot be negative!\n"; - return; - } - - diameter_ = diameter; + vSetDiameter(e,diameter); } double DifShellBase::getDiameter(const Eref& e ) const { - return diameter_; + return vGetDiameter(e); } void DifShellBase::setThickness( const Eref& e, double thickness ) { - if ( thickness < 0.0 ) { - cerr << "Error: DifShell: thickness cannot be negative!\n"; - return; - } - - thickness_ = thickness; + vSetThickness(e,thickness); } double DifShellBase::getThickness(const Eref& e) const { - return thickness_; + return vGetThickness(e); } void DifShellBase::setVolume(const Eref& e, double volume ) { - if ( shapeMode_ != 3 ) - cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n"; - - if ( volume < 0.0 ) { - cerr << "Error: DifShell: volume cannot be negative!\n"; - return; - } - - volume_ = volume; + vSetVolume(e,volume); } double DifShellBase::getVolume(const Eref& e ) const { - return volume_; + return vGetVolume(e); } void DifShellBase::setOuterArea(const Eref& e, double outerArea ) { - if (shapeMode_ != 3 ) - cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n"; - - if ( outerArea < 0.0 ) { - cerr << "Error: DifShell: outerArea cannot be negative!\n"; - return; - } - - outerArea_ = outerArea; + vSetOuterArea(e,outerArea); } double DifShellBase::getOuterArea(const Eref& e ) const { - return outerArea_; + return vGetOuterArea(e); } void DifShellBase::setInnerArea(const Eref& e, double innerArea ) { - if ( shapeMode_ != 3 ) - cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n"; - - if ( innerArea < 0.0 ) { - cerr << "Error: DifShell: innerArea cannot be negative!\n"; - return; - } - - innerArea_ = innerArea; + vSetInnerArea(e,innerArea); } double DifShellBase::getInnerArea(const Eref& e) const { - return innerArea_; + return vGetInnerArea(e); } //////////////////////////////////////////////////////////////////////////////// @@ -416,92 +367,92 @@ void DifShellBase::process( const Eref& e, ProcPtr p ) } void DifShellBase::buffer( - const Eref& e, - double kf, - double kb, - double bFree, - double bBound ) + const Eref& e, + double kf, + double kb, + double bFree, + double bBound ) { vBuffer( e, kf, kb, bFree, bBound ); } void DifShellBase::fluxFromOut(const Eref& e, - double outerC, - double outerThickness ) + double outerC, + double outerThickness ) { vFluxFromOut(e, outerC, outerThickness ); } void DifShellBase::fluxFromIn( - const Eref& e, - double innerC, - double innerThickness ) + const Eref& e, + double innerC, + double innerThickness ) { vFluxFromIn( e, innerC, innerThickness ); } void DifShellBase::influx(const Eref& e, - double I ) + double I ) { vInflux( e, I ); } void DifShellBase::outflux(const Eref& e, - double I ) + double I ) { vOutflux(e, I ); } void DifShellBase::fInflux(const Eref& e, - double I, - double fraction ) + double I, + double fraction ) { vFInflux(e, I, fraction ); } void DifShellBase::fOutflux(const Eref& e, - double I, - double fraction ) + double I, + double fraction ) { vFOutflux(e, I, fraction ); } void DifShellBase::storeInflux(const Eref& e, - double flux ) + double flux ) { vStoreInflux( e, flux ); } void DifShellBase::storeOutflux(const Eref& e, - double flux ) + double flux ) { vStoreOutflux(e, flux ); } void DifShellBase::tauPump(const Eref& e, - double kP, - double Ceq ) + double kP, + double Ceq ) { vTauPump( e, kP, Ceq ); } void DifShellBase::eqTauPump(const Eref& e, - double kP ) + double kP ) { vEqTauPump(e, kP ); } void DifShellBase::mmPump(const Eref& e, - double vMax, - double Kd ) + double vMax, + double Kd ) { vMMPump(e, vMax, Kd ); } void DifShellBase::hillPump(const Eref& e, - double vMax, - double Kd, - unsigned int hill ) + double vMax, + double Kd, + unsigned int hill ) { vHillPump(e, vMax, Kd, hill ); } @@ -509,57 +460,57 @@ void DifShellBase::vSetSolver( const Eref& e, Id hsolve ) {;} void DifShellBase::zombify( Element* orig, const Cinfo* zClass, - Id hsolve ) -{ - if ( orig->cinfo() == zClass ) - return; - unsigned int start = orig->localDataStart(); - unsigned int num = orig->numLocalData(); - if ( num == 0 ) - return; - unsigned int len = 12; - vector< double > data( num * len ); - - unsigned int j = 0; + Id hsolve ) +{ + if ( orig->cinfo() == zClass ) + return; + unsigned int start = orig->localDataStart(); + unsigned int num = orig->numLocalData(); + if ( num == 0 ) + return; + unsigned int len = 12; + vector< double > data( num * len ); + + unsigned int j = 0; - for ( unsigned int i = 0; i < num; ++i ) { - Eref er( orig, i + start ); - const DifShell* ds = - reinterpret_cast< const DifShell* >( er.data() ); - data[j + 0] = ds->getC( er ); - data[j + 1] = ds->getCeq( er ); - data[j + 2] = ds->getD( er ); - data[j + 3] = ds->getValence( er ); - data[j + 4] = ds->getLeak( er ); - data[j + 5] = ds->getShapeMode( er ); - data[j + 6] = ds->getLength( er ); - data[j + 7] = ds->getDiameter( er ); - data[j + 8] = ds->getThickness( er ); - data[j + 9] = ds->getVolume( er ); - data[j + 10] = ds->getOuterArea( er ); - data[j + 11] = ds->getInnerArea( er ); - j += len; - } - orig->zombieSwap( zClass ); - j = 0; - for ( unsigned int i = 0; i < num; ++i ) { - Eref er( orig, i + start ); - const DifShell* ds = - reinterpret_cast< const DifShell* >( er.data() ); - ds->vSetSolver(er,hsolve); - ds->setC(er, data[j+0]); - ds->setCeq(er, data[j + 1]); - ds->setD(er, data[j + 2]); - ds->setValence(er, data[j + 3]); - ds->setLeak(er, data[j + 4]); - ds->setShapeMode(er, data[j + 5]); - ds->setLength(er, data[j + 6]); - ds->setDiameter(er, data[j + 7]); - ds->setThickness(er, data[j + 8]); - ds->setVolume(er, data[j + 9]); - ds->setOuterArea(er, data[j + 10]); - ds->setInnerArea(er, data[j + 11]); - j += len; //?? - } + for ( unsigned int i = 0; i < num; ++i ) { + Eref er( orig, i + start ); + const DifShellBase* ds = + reinterpret_cast< const DifShellBase* >( er.data() ); + data[j + 0] = ds->getC( er ); + data[j + 1] = ds->getCeq( er ); + data[j + 2] = ds->getD( er ); + data[j + 3] = ds->getValence( er ); + data[j + 4] = ds->getLeak( er ); + data[j + 5] = ds->getShapeMode( er ); + data[j + 6] = ds->getLength( er ); + data[j + 7] = ds->getDiameter( er ); + data[j + 8] = ds->getThickness( er ); + data[j + 9] = ds->getVolume( er ); + data[j + 10] = ds->getOuterArea( er ); + data[j + 11] = ds->getInnerArea( er ); + j += len; + } + orig->zombieSwap( zClass ); + j = 0; + for ( unsigned int i = 0; i < num; ++i ) { + Eref er( orig, i + start ); + DifShellBase* ds = + reinterpret_cast< DifShellBase* >( er.data() ); + ds->vSetSolver(er,hsolve); + ds->setC(er, data[j+0]); + ds->setCeq(er, data[j + 1]); + ds->setD(er, data[j + 2]); + ds->setValence(er, data[j + 3]); + ds->setLeak(er, data[j + 4]); + ds->setShapeMode(er, data[j + 5]); + ds->setLength(er, data[j + 6]); + ds->setDiameter(er, data[j + 7]); + ds->setThickness(er, data[j + 8]); + ds->setVolume(er, data[j + 9]); + ds->setOuterArea(er, data[j + 10]); + ds->setInnerArea(er, data[j + 11]); + j += len; //?? + } } diff --git a/moose-core/biophysics/DifShellBase.h b/moose-core/biophysics/DifShellBase.h index 0f413003..db74c536 100644 --- a/moose-core/biophysics/DifShellBase.h +++ b/moose-core/biophysics/DifShellBase.h @@ -16,7 +16,7 @@ class DifShellBase { public: DifShellBase(); - ///////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////// // Dest functions ///////////////////////////////////////////////////////////// void reinit( const Eref & e, ProcPtr p ); @@ -47,7 +47,7 @@ class DifShellBase virtual void vStoreInflux(const Eref& e, double flux ) = 0; virtual void vStoreOutflux(const Eref& e, double flux ) = 0; virtual void vTauPump(const Eref& e, double kP, double Ceq ) = 0; - virtual void rEqTauPump(const Eref& e, double kP ) = 0; + virtual void vEqTauPump(const Eref& e, double kP ) = 0; virtual void vMMPump(const Eref& e, double vMax, double Kd ) = 0; virtual void vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ) = 0; ///////////////////////////////////////////////////////////// @@ -104,12 +104,32 @@ class DifShellBase virtual void vSetLeak(const Eref& e, double leak ) = 0; virtual double vGetLeak(const Eref& e) const = 0; + + virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0; + virtual unsigned int vGetShapeMode(const Eref& e) const = 0; + + virtual void vSetLength(const Eref& e, double length ) = 0; + virtual double vGetLength(const Eref& e) const = 0; + + virtual void vSetDiameter(const Eref& e, double diameter ) = 0; + virtual double vGetDiameter(const Eref& e) const = 0; + + virtual void vSetThickness(const Eref& e, double thickness ) = 0; + virtual double vGetThickness(const Eref& e) const = 0; + + virtual void vSetVolume(const Eref& e, double volume ) = 0; + virtual double vGetVolume(const Eref& e) const = 0; + + virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0; + virtual double vGetOuterArea(const Eref& e) const = 0; + virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0; + virtual double vGetInnerArea(const Eref& e) const = 0; - void vSetSolver( const Eref& e, Id hsolve ); - void zombify( Element* orig, const Cinfo* zClass, - Id hsolve ); + virtual void vSetSolver( const Eref& e, Id hsolve ); + static void zombify( Element* orig, const Cinfo* zClass, + Id hsolve ); static SrcFinfo1< double >* concentrationOut(); static SrcFinfo2< double, double >* innerDifSourceOut(); static SrcFinfo2< double, double >* outerDifSourceOut(); @@ -117,13 +137,6 @@ class DifShellBase private: - unsigned int shapeMode_; - double length_; - double diameter_; - double thickness_; - double volume_; - double outerArea_; - double innerArea_; }; -- GitLab