P& deref_ptr( P *obj) { return *obj; }
+template P& deref_ptr ( P *obj) { return *obj; }
+
+
+class Mesh {
+
+ friend class Cell;
+ friend class Node;
+ friend class FigureEditor;
+
+public:
+ Mesh(void) {
+ // Make sure the reserved value is large enough if a cell is added
+ // in "Divide" when the capacity is insufficient, "cells" might be
+ // relocated including the current Cell (i.e. the value of *this)
+ // calling "Mesh::IncreaseCapacityIfNecessary" (from another
+ // object than Cell, e.g. Mesh) before entering Divide will solve
+ // this issue (solved now).
+ cells.reserve(2);
+ nodes.reserve(500);
+
+ //boundary_polygon = new BoundaryPolygon();
+
+ time = 0.;
+ plugin = 0;
+ };
+ ~Mesh(void) {
+ delete boundary_polygon;
+ /* if (plugin)
+ delete plugin;*/
+ };
+
+ void Clean(void);
+ //void Plane(int xwidth, int ywidth, int nx, int ny, bool randomP=false);
+ Cell &EllipticCell(double xc, double yc, double ra, double rb, int nnodes=10, double rotation=0);
+ Cell &CircularCell(double xc, double yc, double r, int nnodes=10) {
+ return EllipticCell(xc, yc, r, r, nnodes, 0);
+ }
+ Cell &LeafPrimordium(int n, double pet_length);
+ Cell &LeafPrimordium2(int n);
+ Cell *RectangularCell(const Vector ll, const Vector ur, double rotation = 0);
+ void CellFiles(const Vector ll, const Vector ur);
+
+ /* void GMVoutput(ostream &os,
+ const char *codename=0, const char *codever=0,
+ const char *comments=0);*/
+ //void RandPoints(int npoints);
+
+ inline Cell &getCell(int i) {
+ if ((unsigned)i= nodes.size() || i < 0) {
+ // cerr << "Mesh::getNode: Warning. Index " << i << " out of range.\n";
+ // }
+ return *nodes[i];
+ }
+
+ //double Diffusion(void);
+ inline int size(void) {
+ return cells.size();
+ }
+ inline int nnodes(void) {
+ return nodes.size();
+ }
+ //void SortNBLists(void);
+
+ /*template void LoopCells(Op f) {
+ for (vector::iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ f(*i);
+ }
+ }*/
+
+ /*! \brief Calls function f for all Cells f.
+
+ Using this template requires some fiddling with function adaptors bind2nd and mem_fun_ref.
+
+ Example usage for calling a member function on each cell:
+
+ mesh.LoopCells( bind2nd (mem_fun_ref( &Cell::DrawDiffEdges), &canvas ) );
+
+ This calls Cell's member function DrawDiffEdges, taking one
+ argument canvas, on all Cell objects in the current Mesh.
+
+ */
+ template void LoopCells(Op f) {
+ for (vector ::iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ f(**i);
+ }
+ //for_each(cells.begin(),cells.end(),f);
+ }
+
+ template void LoopWalls(Op f) {
+ for (list ::iterator i=walls.begin();
+ i!=walls.end();
+ i++) {
+ f(**i);
+ }
+ }
+
+ // if the amount of cells might increase, during looping, use this template
+ template void LoopCurrentCells(Op f) {
+ vector current_cells = cells;
+ for (vector ::iterator i=current_cells.begin();
+ i!=current_cells.end();
+ i++) {
+ f(**i);
+
+ }
+ //for_each(cells.begin(),cells.end(),f);
+ }
+
+ template void LoopNodes(Op f) {
+ for (vector::iterator i=nodes.begin();
+ i!=nodes.end();
+ i++) {
+ f(**i);
+ }
+ }
+
+ template void RandomlyLoopNodes(Op f) {
+
+ MyUrand r(shuffled_nodes.size());
+ random_shuffle(shuffled_nodes.begin(),shuffled_nodes.end(),r);
+
+ for (vector::const_iterator i=shuffled_nodes.begin();
+ i!=shuffled_nodes.end();
+ i++) {
+ f(*shuffled_nodes[*i]);
+ }
+
+ }
+
+ template void RandomlyLoopCells(Op f) {
+
+ MyUrand r(shuffled_cells.size());
+ random_shuffle(shuffled_cells.begin(),shuffled_cells.end(),r);
+
+ for (vector::const_iterator i=shuffled_cells.begin();
+ i!=shuffled_cells.end();
+ i++) {
+ f(*shuffled_cells[*i]);
+ }
+
+
+ }
+
+ template void LoopCells(Op1 f, Op2 &g) {
+ for (vector::iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ f(**i,g);
+ }
+ }
+
+ template void LoopCells(Op1 f, Op2 &g, Op3 &h) {
+ for (vector::iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ f(**i,g,h);
+ }
+ }
+
+ void DoCellHouseKeeping(void) {
+ vector current_cells = cells;
+ for (vector::iterator i = current_cells.begin();
+ i != current_cells.end();
+ i ++) {
+ plugin->CellHouseKeeping(**i);
+
+ // Call functions of Cell that cannot be called from CellBase, including Division
+ if ((*i)->flag_for_divide) {
+ (*i)->Divide();
+ (*i)->flag_for_divide=false;
+ }
+ }
+ }
+/* template void ExtractFromCells(Op1 f, Cont res) {
+ for (vector::iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ *(res++) = ( f(*i) );
+ }
+ }*/
+
+ // Apply "f" to cell i
+ // i.e. this is an adapter which allows you to call a function
+ // operating on Cell on its numeric index index
+ template void cell_index_adapter(Op f,int i) {
+ f(cells[i]);
+ }
+
+ double DisplaceNodes(void);
+
+ void BoundingBox(Vector &LowerLeft, Vector &UpperRight);
+ int NEqs(void) { int nwalls = walls.size();
+ int ncells =cells.size();
+ int nchems = Cell::NChem();
+
+ // two eqs per chemical for each walls, and one eq per chemical for each cell
+ // This is for generality. For a specific model you may optimize
+ // this by removing superfluous (empty) equations.
+ int neqs = 2 * nwalls * nchems + ncells * nchems;
+
+ return neqs;
+ }
+ void IncreaseCellCapacityIfNecessary(void) {
+
+ return;
+ // cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";
+ // make sure we always have enough space
+ // to have each cell divide at least once
+ //
+ // Note that we must do this, because Cell::Divide pushes a new Cell
+ // onto Mesh::cells. As a result, Mesh::cells might be relocated
+ // if we are _within_ a Cell object: i.e. pointer "this" will be changed!!
+ //
+ // An alternative solution could be to make "Mesh::cells" a list,
+ // but this won't work because we need random access for
+ // the Monte Carlo algorithm.
+
+ if (2*cells.size()>cells.capacity()) {
+ cerr << "Increasing capacity to " << 2*cells.capacity() << endl;
+ cerr << "Current capacity is " << cells.capacity() << endl;
+ cells.reserve(cells.capacity()*2);
+ }
+ }
+
+ void ReserveMoreCells(int n) {
+ if (nodes.size()+n>nodes.capacity()) {
+ nodes.reserve(size()+n);
+ }
+ }
+ double Area(void);
+ double MeanArea(void) {
+ double sum=0.;
+ for (vector::const_iterator i=cells.begin();
+ i!=cells.end();
+ i++) {
+ sum+=(*i)->Area();
+ }
+ return sum/(double)NCells();
+ }
+
+ void SetBaseArea(void);
+ int NCells(void) const {
+ return cells.size();
+ }
+ inline int NNodes(void) const {
+ return nodes.size();
+ }
+ void PrintQueue(ostream &os) {
+ while (!node_insertion_queue.empty()) {
+ os << node_insertion_queue.front() << endl;
+ node_insertion_queue.pop();
+ }
+ //copy (node_insertion_queue.begin(),node_insertion_queue.end(),ostream_iterator(cerr, " "));
+ }
+
+ void InsertNodes(void) {
+ // insert the nodes in the insertion queue
+ while (!node_insertion_queue.empty()) {
+
+ //cerr << node_insertion_queue.front() << endl;
+ InsertNode(node_insertion_queue.front());
+ node_insertion_queue.pop();
+ }
+
+ }
+
+ void Clear();
+
+ // template ReactDiffuse(const double D, ReactFunction& react) {
+ void ReactDiffuse( double delta_t = 1 );
+ //void Diffuse(const double D);
+ double SumChemical(int ch);
+ void SetChemical(int ch, double value) {
+ for (vector::iterator c=cells.begin();
+ c!=cells.end();
+ c++) {
+ (*c)->chem[ch]=value;
+ }
+ }
+
+ // used for interacing with ODE-solvers (e.g. NRCRungeKutta)
+ void setValues(double x, double *y);
+ double *getValues(int *neqs);
+ void Derivatives(double *derivs);
+#ifdef QTGRAPHICS
+ inline void DrawBoundary(QGraphicsScene *c) {
+ boundary_polygon->Draw(c);
+ }
+ void DrawNodes(QGraphicsScene *c) const;
+
+#endif
+ double max_chem;
+
+ void XMLSave(const char *docname, xmlNode *settings=0) const;
+ void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);
+ void XMLReadPars(const xmlNode * root_node);
+ void XMLReadGeometry(const xmlNode *root_node);
+ void XMLReadSimtime(const xmlNode *root_node);
+ void XMLReadNodes(xmlNode *cur);
+ void XMLReadCells(xmlNode *cur);
+ void XMLParseTree(const xmlNode * root_node);
+ void XMLReadWalls(xmlNode *cur, vector *tmp_cells);
+ void XMLReadWallsToCells(xmlNode *root, vector *tmp_walls);
+ void XMLReadNodeSets(xmlNode *root);
+ void XMLReadNodeSetsToNodes(xmlNode *root);
+ void PerturbChem(int chemnum, double range);
+ void CleanUpCellNodeLists(void);
+ void CleanUpWalls(void);
+ void CutAwayBelowLine( Vector startpoint, Vector endpoint );
+ void CutAwaySAM(void);
+ void RepairBoundaryPolygon(void);
+ void Rotate(double angle, Vector center);
+ void PrintWallList( void );
+ void TestIllegalWalls(void);
+ Vector FirstConcMoment(int chem);
+ inline Vector Centroid(void) {
+ return boundary_polygon->Centroid();
+ }
+
+ inline Vector Offset(void) {
+ return boundary_polygon->Offset();
+ }
+
+ inline double Factor(void) {
+ return boundary_polygon->Factor();
+ }
+
+ void DeleteLooseWalls(void);
+ void FitLeafToCanvas(double width, double height);
+ void AddNodeSet(NodeSet *node_set) {
+ node_sets.push_back(node_set);
+ }
+
+ void CleanChemicals(const vector &clean_chem, const vector &clean_transporters);
+ void RandomizeChemicals(const vector &max_chem, const vector &max_transporters);
+ inline double getTime(void) const { return time; }
+ string getTimeHours(void) const;
+ inline void setTime(double t) { time = t; }
+ double CalcProtCellsWalls(int ch) const;
+ void SettoInitVals(void);
+ QVector VertexAngles(void);
+ QVector< QPair > VertexAnglesValues(void);
+ void SetSimPlugin(SimPluginInterface *new_plugin) {
+ /* if (plugin)
+ delete plugin;*/
+ plugin=new_plugin;
+ }
+ QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
+ void StandardInit(void);
+
+ Node* findNextBoundaryNode(Node*);
+
+private:
+
+ // Data members
+ vector cells;
+ vector nodes;
+ list walls; // we need to erase elements from this container frequently, hence a list.
+public:
+ vector node_sets;
+private:
+ vector shuffled_nodes;
+ vector| shuffled_cells;
+ unique_queue node_insertion_queue;
+ BoundaryPolygon *boundary_polygon;
+ double time;
+ SimPluginInterface *plugin;
+
+ // Private member functions
+ void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
+ void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list::iterator ins_pos);
+ void InsertNode(Edge &e);
+ inline Node *AddNode(Node *n) {
+ nodes.push_back(n);
+ shuffled_nodes.push_back(n);
+ n->m=this;
+ return n;
+ }
+
+ inline Cell *AddCell(Cell *c) {
+ cells.push_back(c);
+ shuffled_cells.push_back(c);
+ //cerr << "Shuffled cell indices: ";
+ /*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator(cerr," "));
+ cerr << endl;*/
+ c->m=this;
+ return c;
+ }
+
+ //int Delaunay(void);
+ void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
+ double *xc,double *yc,double *r);
+
+
+ // void RenumberCells(void);
+
+};
+#endif
| | | | | | | | | | | | | |