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) {
if ((*i)->division_axis) {
(*i)->DivideOverAxis(*(*i)->division_axis);
delete (*i)->division_axis;
(*i)->division_axis = 0;
} else {
(*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);
void CleanTransporters(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
| | | | | | | | | | | | | |