VPTissue Reference Manual
Cell.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2011-2016 Universiteit Antwerpen
3  *
4  * Licensed under the EUPL, Version 1.1 or as soon they will be approved by
5  * the European Commission - subsequent versions of the EUPL (the "Licence");
6  * You may not use this work except in compliance with the Licence.
7  * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the Licence is distributed on an "AS IS" basis,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the Licence for the specific language governing
13  * permissions and limitations under the Licence.
14  */
20 #include "Cell.h"
21 
22 #include "AreaMoment.h"
23 #include "CellAttributes.h"
24 #include "Edge.h"
25 #include "GeoData.h"
26 #include "Node.h"
27 
28 #include "math/array3.h"
29 #include "math/constants.h"
30 #include "util/container/CircularIterator.h"
31 #include "util/container/ConstCircularIterator.h"
32 
33 #include <iomanip>
34 #include <list>
35 #include <string>
36 #include <tuple>
37 #include <vector>
38 
39 namespace SimPT_Sim { class Node; }
40 namespace SimPT_Sim { class Wall; }
41 
42 using namespace std;
43 using namespace boost::property_tree;
44 using boost::optional;
45 using namespace SimPT_Sim::Container;
46 using namespace SimPT_Sim::Util;
47 using namespace SimPT_Sim;
48 
49 namespace {
50 
56 bool is_intersecting(const array<double, 3>& v1, const array<double, 3>& v2,
57  const array<double, 3>& v3, const array<double, 3>& v4)
58 {
59  double const denominator = (v4[1] - v3[1]) * (v2[0] - v1[0]) - (v4[0] - v3[0]) * (v2[1] - v1[1]);
60  double const ua = ((v4[0] - v3[0]) * (v1[1] - v3[1]) - (v4[1] - v3[1]) * (v1[0] - v3[0]))
61  / denominator;
62  double const ub = ((v2[0] - v1[0]) * (v1[1] - v3[1]) - (v2[1] - v1[1]) * (v1[0] - v3[0]))
63  / denominator;
64 
65  return (0.0 < ua) && (ua < 1.0) && (0.0 < ub) && (ub < 1.0);
66 }
67 
68 } // namespace anonymous
69 
70 namespace SimPT_Sim {
71 
72 Cell::Cell(int index, unsigned int chem_count)
73  : CellAttributes(chem_count), m_index(index),
74  m_dirty_area(true), m_dirty_centroid(true), m_dirty_geo_data(true)
75 {
76 }
77 
78 void Cell::AddWall(Wall* w)
79 {
80  // Add Wall to Cell's list (already registered with Mesh by Mesh::BuildWall)
81  m_walls.push_back(w);
82 }
83 
84 void Cell::CalculateArea() const
85 {
86  double a = 0.0;
87 
88  auto cit = make_const_circular(m_nodes, m_nodes.begin());
89  do {
90  const double i_x = (*(*cit))[0];
91  const double i_y = (*(*cit))[1];
92  const double n_x = (*(*next(cit)))[0];
93  const double n_y = (*(*next(cit)))[1];
94 
95  a += i_x * n_y - n_x * i_y;
96  } while (++cit != m_nodes.begin());
97 
98  m_geo_data.m_area = abs(a) / 2.0;
99  m_dirty_area = false;
100 }
101 
102 void Cell::CalculateCentroid() const
103 {
104  double a = 0.0; // area
105  double x = 0.0; // centroid - x component
106  double y = 0.0; // centroid - y component
107 
108  auto cit = make_const_circular(m_nodes, m_nodes.begin());
109  do {
110  const double i_x = (*(*cit))[0];
111  const double i_y = (*(*cit))[1];
112  const double n_x = (*(*next(cit)))[0];
113  const double n_y = (*(*next(cit)))[1];
114 
115  a += i_x * n_y - n_x * i_y;
116  x += (n_x + i_x) * (i_x * n_y - n_x * i_y);
117  y += (n_y + i_y) * (i_x * n_y - n_x * i_y);
118 
119  } while (++cit != m_nodes.begin());
120 
121  m_geo_data.m_area = 0.5 * abs(a);
122  m_geo_data.m_centroid = array<double, 3> {{ x / (3.0 * a), y / (3.0 * a), 0.0}};
123  m_dirty_area = false;
124  m_dirty_centroid = false;
125 }
126 
127 void Cell::CalculateGeoData() const
128 {
129  double a = 0.0; // area (actually twice the area)
130  double x = 0.0; // centroid - x component
131  double y = 0.0; // centroid - y component
132  double xx = 0.0; // area moment of inertia - xx component
133  double xy = 0.0; // area moment of inertia - xy component
134  double yy = 0.0; // area moment of inertia - yy component
135 
136  auto cit = make_const_circular(m_nodes, m_nodes.begin());
137  do {
138  const double i_x = (*(*cit))[0];
139  const double i_y = (*(*cit))[1];
140  const double n_x = (*(*next(cit)))[0];
141  const double n_y = (*(*next(cit)))[1];
142 
143  a += i_x * n_y - n_x * i_y;
144  x += (n_x + i_x) * (i_x * n_y - n_x * i_y);
145  y += (n_y + i_y) * (i_x * n_y - n_x * i_y);
146  xx += (i_x * i_x + n_x * i_x + n_x * n_x) * (i_x * n_y - n_x * i_y);
147  xy += (n_x * i_y - i_x * n_y) * (i_x * (2 * i_y + n_y) + n_x * (i_y + 2 * n_y));
148  yy += (i_x * n_y - n_x * i_y) * (i_y * i_y + n_y * i_y + n_y * n_y);
149 
150  } while (++cit != m_nodes.begin());
151 
152  m_geo_data.m_area = abs(a) / 2.0;
153  m_geo_data.m_centroid = array<double, 3> {{ x / (3.0 * a), y / (3.0 * a), 0.0}};
154  m_geo_data.m_area_moment = AreaMoment(xx / 12.0, xy / 24.0, yy / 12.0);
155 
156  m_dirty_area = false;
157  m_dirty_centroid = false;
158  m_dirty_geo_data = false;
159 }
160 
161 ptree Cell::GeometryToPtree() const
162 {
163  ptree ret;
164  ret.put("id", m_index);
165  ret.put("area", GetArea());
166  ret.put("at_boundary", HasBoundaryWall());
167 
168  for (const auto& node : m_nodes) {
169  ret.add("node_array.node", node->GetIndex());
170  }
171  for (const auto& wall : m_walls) {
172  ret.add("wall_array.wall", wall->GetIndex());
173  }
174 
175  return ret;
176 }
177 
178 double Cell::GetArea() const
179 {
180  if (m_dirty_area) {
181  CalculateArea();
182  }
183  return m_geo_data.m_area;
184 }
185 
186 double Cell::GetSignedArea() const
187 {
188  // Should be handled similar to area but first we need to test if it solves
189  // the problem of cell inversion; the sole reason it's been introduced.
190  double a = 0.0;
191 
192  auto cit = make_const_circular(m_nodes, m_nodes.begin());
193  do {
194  const double i_x = (*(*cit))[0];
195  const double i_y = (*(*cit))[1];
196  const double n_x = (*(*next(cit)))[0];
197  const double n_y = (*(*next(cit)))[1];
198 
199  a += i_x * n_y - n_x * i_y;
200  } while (++cit != m_nodes.begin());
201 
202  return 0.5 * a;
203 }
204 
205 array<double, 3> Cell::GetCentroid() const
206 {
207  if (m_dirty_centroid) {
208  CalculateCentroid();
209  }
210  return m_geo_data.m_centroid;
211 }
212 
213 double Cell::GetCircumference() const
214 {
215  double circumference = 0.0;
216  auto cit = make_const_circular(m_nodes);
217  do {
218  const double dx = (*(*next(cit)))[0] - (*(*cit))[0];
219  const double dy = (*(*next(cit)))[1] - (*(*cit))[1];
220  circumference += sqrt(dx * dx + dy * dy);
221  } while (++cit != m_nodes.begin());
222  return circumference;
223 }
224 
225 GeoData Cell::GetGeoData() const
226 {
227  if (m_dirty_geo_data) {
228  CalculateGeoData();
229  }
230  return m_geo_data;
231 }
232 
233 double Cell::GetSumTransporters(unsigned int ch) const
234 {
235  double sum = 0.0;
236  for (const auto& wall : m_walls) {
237  assert( (wall->GetC1() == this || wall->GetC2() == this) && "Error in cell <-> wall!" );
238  sum += (wall->GetC1() == this) ? wall->GetTransporters1(ch) : wall->GetTransporters2(ch);
239  }
240  return sum;
241 }
242 
243 bool Cell::HasEdge(const Edge& edge) const
244 {
245  bool status = false;
246  auto cit = make_const_circular(m_nodes, m_nodes.begin());
247  do {
248  if (Edge(*cit, *next(cit)) == edge) {
249  status = true;
250  break;
251  }
252  } while(++cit != m_nodes.begin());
253  return status;
254 }
255 
256 bool Cell::HasBoundaryWall() const
257 {
258  bool status = false;
259  for (const auto& wall : m_walls) {
260  if (wall->IsAtBoundary()) {
261  status = true;
262  break;
263  }
264  }
265  return status;
266 }
267 
268 bool Cell::HasNeighborOfTypeZero() const
269 {
270  int prod = 1;
271  for (const auto& wall : m_walls) {
272 
273 /* prod *= ( wall->GetC1()->GetIndex() != this->GetIndex()
274  ? ( wall->GetC1()->GetCellType() ) : ( wall->GetC2()->GetCellType() ) );
275 
276 */
277  if ( wall->GetC1()->GetIndex() != this->GetIndex() ) {
278  if ( wall->GetC1()->GetIndex() != -1 ) {
279  prod *= wall->GetC1()->GetCellType();
280  }
281  } else {
282  if ( wall->GetC2()->GetIndex() != -1 ) {
283  prod *= wall->GetC2()->GetCellType();
284  }
285  }
286  }
287  return ( prod == 0 );
288 }
289 
290 bool Cell::IsWallNeighbor(Cell* cell) const
291 {
292  bool status = false;
293  for (const auto& wall : m_walls) {
294  if ((cell != this) && (wall->GetC1() == cell || wall->GetC2() == cell)) {
295  status = true;
296  break;
297  }
298  }
299  return status;
300 }
301 
302 void Cell::Move(const array<double, 3>& a)
303 {
304  for (auto& node : m_nodes) {
305  *node += a;
306  }
307 }
308 
309 bool Cell::MoveSelfIntersects(Node* moving_node, array<double, 3> new_pos)
310 {
311  const auto moving_node_it = find(m_nodes.begin(), m_nodes.end(), moving_node);
312 
313  // Only makes sense if there are nodes and moving_node is among them
314  if (!m_nodes.empty() && moving_node_it != m_nodes.end()) {
315 
316  array<double, 3> neighbor_of_moving[2];
317  const auto moving_it = make_const_circular(m_nodes, moving_node_it);
318  neighbor_of_moving[0] = *(*next(moving_it));
319  neighbor_of_moving[1] = *(*prev(moving_it));
320 
321  // Compare the two new edges against edge starting
322  // at node *it, with it looping over m_nodes
323  auto it = make_const_circular(m_nodes);
324  do {
325  // loop over the two neighbors of moving node
326  for (int j = 0; j < 2; j++) {
327  // do not compare to self
328  if (*it == moving_node || *next(it) == moving_node) {
329  continue;
330  }
331  array<double, 3> v3 = *(*it);
332  array<double, 3> v4 = *(*next(it));
333  if (is_intersecting(new_pos, neighbor_of_moving[j], v3, v4)) {
334  return true;
335  }
336  }
337  } while(++it != m_nodes.begin());
338  }
339  return false;
340 }
341 
342 ostream& Cell::Print(ostream& os) const
343 {
344  os << "Cell: [ index = " << m_index << " ]" << endl;
345 
346  os << "Nodes: { " << endl << "\t ";
347  for (auto const& node : m_nodes) {
348  os << node->GetIndex() << " && ";
349  }
350  os << endl <<"} " << endl;
351  os << "Walls: { " << endl << "\t ";
352  for (auto const& wall : m_walls) {
353  wall->Print(os);
354  os << ", ";
355  }
356  os << endl << "} " << endl;
357 
358  CellAttributes::Print(os);
359 
360  return os;
361 }
362 
363 void Cell::ReadPtree(const ptree& cell_pt)
364 {
365  if ( !IsBoundaryPolygon() ) {
366  const ptree& attributes_pt = cell_pt.get_child("attributes");
367  CellAttributes::ReadPtree(attributes_pt);
368  }
369 }
370 
371 void Cell::ReassignWall(Wall* w, Cell* to)
372 {
373  // Remove wall from this cell's list
374  m_walls.remove(w);
375  // Reassign it to the "to" cell
376  to->m_walls.push_back(w);
377 }
378 
379 void Cell::SetGeoDirty()
380 {
381  m_dirty_area = true;
382  m_dirty_centroid = true;
383  m_dirty_geo_data = true;
384 }
385 
386 void Cell::SetTransporters(unsigned int ch, double conc)
387 {
388  for (auto const& wall : m_walls) {
389  assert( (wall->GetC1() == this || wall->GetC2() == this) && "Error in cell <-> wall!" );
390 
391  if (wall->GetC1() == this) {
392  wall->SetTransporters1(ch,conc);
393  } else {
394  wall->SetTransporters2(ch, conc);
395  }
396  }
397 }
398 
399 void Cell::SetTransporters(unsigned int ch, double , double lat) //DDV: for now this makes sense only if conc = 1
400 {
401  for (auto w = GetWalls().begin(); w != GetWalls().end(); ++w) {
402 
403  std::array<double, 3> ref{{1., 0., 0.}};
404  std::array<double, 3> wa = *((*w)->GetN2()) - *((*w)->GetN1());
405  double t2 = SignedAngle(wa, ref);
406 
407  if (this->GetCellType() == 1)
408  {
409  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
410  { //UPPER WALL
411  (*w)->SetTransporters1(ch, 1.);
412  (*w)->SetTransporters2(ch, 0.);
413  }
414  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
415  { //LOWER WALL
416  (*w)->SetTransporters1(ch, 0.);
417  (*w)->SetTransporters2(ch, 1.);
418  }
419  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
420  {
421  if ((*w)->GetC1() == this)
422  { //LEFT WALL WRT C1
423  (*w)->SetTransporters1(ch, 0.);
424  (*w)->SetTransporters2(ch, 0.);
425  }
426  else
427  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
428  (*w)->SetTransporters1(ch, 0.);
429  (*w)->SetTransporters2(ch, lat);
430  }
431  }
432  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
433  {
434  if ((*w)->GetC1() == this)
435  { //RIGHT WALL WRT C1
436  (*w)->SetTransporters1(ch, lat);
437  (*w)->SetTransporters2(ch, 0.);
438  }
439  else
440  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
441  (*w)->SetTransporters1(ch, 0.);
442  (*w)->SetTransporters2(ch, 0.);
443  }
444  }
445  } else if (this->GetCellType() == 2) {
446 
447  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
448  { //UPPER WALL
449  (*w)->SetTransporters1(ch, 1.);
450  (*w)->SetTransporters2(ch, 0.);
451  }
452  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
453  { //LOWER WALL
454  (*w)->SetTransporters1(ch, 0.);
455  (*w)->SetTransporters2(ch, 1.);
456  }
457  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
458  {
459  if ((*w)->GetC1() == this)
460  { //LEFT WALL WRT C1
461  (*w)->SetTransporters1(ch, 0.);
462  (*w)->SetTransporters2(ch, lat);
463  }
464  else
465  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
466  (*w)->SetTransporters1(ch, 0.);
467  (*w)->SetTransporters2(ch, lat);
468  }
469  }
470  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
471  {
472  if ((*w)->GetC1() == this)
473  { //RIGHT WALL WRT C1
474  (*w)->SetTransporters1(ch, lat);
475  (*w)->SetTransporters2(ch, 0.);
476  }
477  else
478  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
479  (*w)->SetTransporters1(ch, lat);
480  (*w)->SetTransporters2(ch, 0.);
481  }
482  }
483  } else if (this->GetCellType() == 3) {
484 
485  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
486  { //UPPER WALL
487  (*w)->SetTransporters1(ch, 0.);
488  (*w)->SetTransporters2(ch, 1.);
489  }
490  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
491  { //LOWER WALL
492  if (this->HasNeighborOfTypeZero())
493  {
494  (*w)->SetTransporters1(ch, 1.);
495  (*w)->SetTransporters2(ch, 1.);
496  }
497  else
498  {
499  (*w)->SetTransporters1(ch, 1.);
500  (*w)->SetTransporters2(ch, 0.);
501  }
502  }
503  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
504  {
505  if ((*w)->GetC1() == this)
506  { //LEFT WALL WRT C1
507  (*w)->SetTransporters1(ch, 0.);
508  (*w)->SetTransporters2(ch, lat);
509  }
510  else
511  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
512  (*w)->SetTransporters1(ch, 0.);
513  (*w)->SetTransporters2(ch, lat);
514  }
515  }
516  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
517  {
518  if ((*w)->GetC1() == this)
519  { //RIGHT WALL WRT C1
520  (*w)->SetTransporters1(ch, lat);
521  (*w)->SetTransporters2(ch, 0.);
522  }
523  else
524  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
525  (*w)->SetTransporters1(ch, lat);
526  (*w)->SetTransporters2(ch, 0.);
527  }
528  }
529  } else if (this->GetCellType() == 4) {
530 
531  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
532  { //UPPER WALL
533  (*w)->SetTransporters1(ch, 0.);
534  (*w)->SetTransporters2(ch, 1.);
535  }
536  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
537  { //LOWER WALL
538  if (this->HasNeighborOfTypeZero())
539  {
540  (*w)->SetTransporters1(ch, 1.);
541  (*w)->SetTransporters2(ch, 1.);
542  }
543  else
544  {
545  (*w)->SetTransporters1(ch, 1.);
546  (*w)->SetTransporters2(ch, 0.);
547  }
548  }
549  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
550  {
551  if ((*w)->GetC1() == this)
552  { //LEFT WALL WRT C1
553  (*w)->SetTransporters1(ch, 0.);
554  (*w)->SetTransporters2(ch, lat);
555  }
556  else
557  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
558  (*w)->SetTransporters1(ch, 0.);
559  (*w)->SetTransporters2(ch, 0.);
560  }
561  }
562  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
563  {
564  if ((*w)->GetC1() == this)
565  { //RIGHT WALL WRT C1
566  (*w)->SetTransporters1(ch, 0.);
567  (*w)->SetTransporters2(ch, 0.);
568  }
569  else
570  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
571  (*w)->SetTransporters1(ch, lat);
572  (*w)->SetTransporters2(ch, 0.);
573  }
574  }
575  } else if (this->GetCellType() == 5) {
576 
577  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
578  { //UPPER WALL
579  (*w)->SetTransporters1(ch, 0.);
580  (*w)->SetTransporters2(ch, 1.);
581  }
582  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
583  { //LOWER WALL
584  if (this->HasNeighborOfTypeZero())
585  {
586  (*w)->SetTransporters1(ch, 1.);
587  (*w)->SetTransporters2(ch, 1.);
588  }
589  else
590  {
591  (*w)->SetTransporters1(ch, 1.);
592  (*w)->SetTransporters2(ch, 0.);
593  }
594  }
595  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
596  {
597  if ((*w)->GetC1() == this)
598  { //LEFT WALL WRT C1
599  (*w)->SetTransporters1(ch, 0.);
600  (*w)->SetTransporters2(ch, 0.);
601  }
602  else
603  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
604  (*w)->SetTransporters1(ch, 0.);
605  (*w)->SetTransporters2(ch, 0.);
606  }
607  }
608  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
609  {
610  if ((*w)->GetC1() == this)
611  { //RIGHT WALL WRT C1
612  (*w)->SetTransporters1(ch, 0.);
613  (*w)->SetTransporters2(ch, 0.);
614  }
615  else
616  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
617  (*w)->SetTransporters1(ch, 0.);
618  (*w)->SetTransporters2(ch, 0.);
619  }
620  }
621  } else if (this->GetCellType() == 6) {
622  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
623  { //UPPER WALL
624  (*w)->SetTransporters1(ch, 0.);
625  (*w)->SetTransporters2(ch, 1.);
626  }
627  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
628  { //LOWER WALL
629  if (this->HasNeighborOfTypeZero())
630  {
631  (*w)->SetTransporters1(ch, 1.);
632  (*w)->SetTransporters2(ch, 1.);
633  }
634  else
635  {
636  (*w)->SetTransporters1(ch, 1.);
637  (*w)->SetTransporters2(ch, 0.);
638  }
639  }
640  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
641  {
642  if ((*w)->GetC1() == this)
643  { //LEFT WALL WRT C1
644  (*w)->SetTransporters1(ch, 0.);
645  (*w)->SetTransporters2(ch, 0.);
646  }
647  else
648  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
649  (*w)->SetTransporters1(ch, lat);
650  (*w)->SetTransporters2(ch, 0.);
651  }
652  }
653  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
654  {
655  if ((*w)->GetC1() == this)
656  { //RIGHT WALL WRT C1
657  (*w)->SetTransporters1(ch, 0.);
658  (*w)->SetTransporters2(ch, lat);
659  }
660  else
661  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
662  (*w)->SetTransporters1(ch, 0.);
663  (*w)->SetTransporters2(ch, 0.);
664  }
665  }
666  } else if (this->GetCellType() == 7) {
667 
668  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
669  { //UPPER WALL
670  (*w)->SetTransporters1(ch, 0.);
671  (*w)->SetTransporters2(ch, 1.);
672  }
673  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
674  { //LOWER WALL
675  if (this->HasNeighborOfTypeZero())
676  {
677  (*w)->SetTransporters1(ch, 1.);
678  (*w)->SetTransporters2(ch, 1.);
679  }
680  else
681  {
682  (*w)->SetTransporters1(ch, 1.);
683  (*w)->SetTransporters2(ch, 0.);
684  }
685  }
686  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
687  {
688  if ((*w)->GetC1() == this)
689  { //LEFT WALL WRT C1
690  (*w)->SetTransporters1(ch, lat);
691  (*w)->SetTransporters2(ch, 0.);
692  }
693  else
694  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
695  (*w)->SetTransporters1(ch, lat);
696  (*w)->SetTransporters2(ch, 0.);
697  }
698  }
699  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
700  {
701  if ((*w)->GetC1() == this)
702  { //RIGHT WALL WRT C1
703  (*w)->SetTransporters1(ch, 0.);
704  (*w)->SetTransporters2(ch, lat);
705  }
706  else
707  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
708  (*w)->SetTransporters1(ch, 0.);
709  (*w)->SetTransporters2(ch, lat);
710  }
711  }
712  } else if (this->GetCellType() == 8) {
713 
714  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
715  { //UPPER WALL
716  (*w)->SetTransporters1(ch, 1.);
717  (*w)->SetTransporters2(ch, 0.);
718  }
719  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
720  { //LOWER WALL
721  (*w)->SetTransporters1(ch, 0.);
722  (*w)->SetTransporters2(ch, 1.);
723  }
724  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
725  {
726  if ((*w)->GetC1() == this)
727  { //LEFT WALL WRT C1
728  (*w)->SetTransporters1(ch, lat);
729  (*w)->SetTransporters2(ch, 0.);
730  }
731  else
732  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
733  (*w)->SetTransporters1(ch, lat);
734  (*w)->SetTransporters2(ch, 0.);
735  }
736  }
737  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
738  {
739  if ((*w)->GetC1() == this)
740  { //RIGHT WALL WRT C1
741  (*w)->SetTransporters1(ch, 0.);
742  (*w)->SetTransporters2(ch, lat);
743  }
744  else
745  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
746  (*w)->SetTransporters1(ch, 0.);
747  (*w)->SetTransporters2(ch, lat);
748  }
749  }
750  } else if (this->GetCellType() == 9) {
751 
752  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
753  { //UPPER WALL
754  (*w)->SetTransporters1(ch, 1.);
755  (*w)->SetTransporters2(ch, 0.);
756  }
757  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
758  { //LOWER WALL
759  (*w)->SetTransporters1(ch, 0.);
760  (*w)->SetTransporters2(ch, 1.);
761  }
762 
763  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
764  {
765  if ((*w)->GetC1() == this)
766  { //LEFT WALL WRT C1
767  (*w)->SetTransporters1(ch, lat);
768  (*w)->SetTransporters2(ch, 0.);
769  }
770  else
771  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
772  (*w)->SetTransporters1(ch, 0.);
773  (*w)->SetTransporters2(ch, 0.);
774  }
775  }
776  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
777  {
778  if ((*w)->GetC1() == this)
779  { //RIGHT WALL WRT C1
780  (*w)->SetTransporters1(ch, 0.);
781  (*w)->SetTransporters2(ch, 0.);
782  }
783  else
784  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
785  (*w)->SetTransporters1(ch, 0.);
786  (*w)->SetTransporters2(ch, lat);
787  }
788  }
789  }
790  }
791 }
792 
793 void Cell::SetTransportersLight(unsigned int ch, double , double lat) //DDV: for now this makes sense only if conc = 1
794 {
795  for (list<Wall *>::iterator w = GetWalls().begin(); w != GetWalls().end(); w++) {
796 
797  std::array<double, 3> ref{{1., 0., 0.}};
798  std::array<double, 3> wa = *((*w)->GetN2()) - *((*w)->GetN1());
799  double t2 = SignedAngle(wa, ref);
800 
801  if (this->GetCellType() == 1)
802  {
803  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
804  { //UPPER WALL
805  (*w)->SetTransporters1(ch, 1.);
806  (*w)->SetTransporters2(ch, 0.);
807  }
808  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
809  { //LOWER WALL
810  (*w)->SetTransporters1(ch, 0.);
811  (*w)->SetTransporters2(ch, 1.);
812  }
813  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
814  {
815  if ((*w)->GetC1() == this)
816  { //LEFT WALL WRT C1
817  (*w)->SetTransporters1(ch, 0.);
818  (*w)->SetTransporters2(ch, 0.);
819  }
820  else
821  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
822  (*w)->SetTransporters1(ch, 0.);
823  (*w)->SetTransporters2(ch, lat);
824  }
825  }
826  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
827  {
828  if ((*w)->GetC1() == this)
829  { //RIGHT WALL WRT C1
830  (*w)->SetTransporters1(ch, lat);
831  (*w)->SetTransporters2(ch, 0.);
832  }
833  else
834  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
835  (*w)->SetTransporters1(ch, 0.);
836  (*w)->SetTransporters2(ch, 0.);
837  }
838  }
839  }
840  else if (this->GetCellType() == 2) {
841 
842  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
843  { //UPPER WALL
844  (*w)->SetTransporters1(ch, 0.);
845  (*w)->SetTransporters2(ch, 1.);
846  }
847  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
848  { //LOWER WALL
849  if (this->HasNeighborOfTypeZero())
850  {
851  (*w)->SetTransporters1(ch, 1.);
852  (*w)->SetTransporters2(ch, 1.);
853  }
854  else
855  {
856  (*w)->SetTransporters1(ch, 1.);
857  (*w)->SetTransporters2(ch, 0.);
858  }
859  }
860  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
861  {
862  if ((*w)->GetC1() == this)
863  { //LEFT WALL WRT C1
864  (*w)->SetTransporters1(ch, 0.);
865  (*w)->SetTransporters2(ch, lat);
866  }
867  else
868  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
869  (*w)->SetTransporters1(ch, 0.);
870  (*w)->SetTransporters2(ch, 0.);
871  }
872  }
873  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
874  {
875  if ((*w)->GetC1() == this)
876  { //RIGHT WALL WRT C1
877  (*w)->SetTransporters1(ch, 0.);
878  (*w)->SetTransporters2(ch, 0.);
879  }
880  else
881  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
882  (*w)->SetTransporters1(ch, lat);
883  (*w)->SetTransporters2(ch, 0.);
884  }
885  }
886  }
887  else if (this->GetCellType() == 3) {
888  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
889  { //UPPER WALL
890  (*w)->SetTransporters1(ch, 0.);
891  (*w)->SetTransporters2(ch, 1.);
892  }
893  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
894  { //LOWER WALL
895  if (this->HasNeighborOfTypeZero())
896  {
897  (*w)->SetTransporters1(ch, 1.);
898  (*w)->SetTransporters2(ch, 1.);
899  }
900  else
901  {
902  (*w)->SetTransporters1(ch, 1.);
903  (*w)->SetTransporters2(ch, 0.);
904  }
905  }
906  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
907  {
908  if ((*w)->GetC1() == this)
909  { //LEFT WALL WRT C1
910  (*w)->SetTransporters1(ch, 0.);
911  (*w)->SetTransporters2(ch, 0.);
912  }
913  else
914  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
915  (*w)->SetTransporters1(ch, lat);
916  (*w)->SetTransporters2(ch, 0.);
917  }
918  }
919  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
920  {
921  if ((*w)->GetC1() == this)
922  { //RIGHT WALL WRT C1
923  (*w)->SetTransporters1(ch, 0.);
924  (*w)->SetTransporters2(ch, lat);
925  }
926  else
927  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
928  (*w)->SetTransporters1(ch, 0.);
929  (*w)->SetTransporters2(ch, 0.);
930  }
931  }
932  else if (this->GetCellType() == 4) {
933 
934  if ((t2 > (-pi() / 4)) && (t2 <= (pi() / 4)))
935  { //UPPER WALL
936  (*w)->SetTransporters1(ch, 1.);
937  (*w)->SetTransporters2(ch, 0.);
938  }
939  else if ((t2 <= (-pi() * 3 / 4)) || (t2 > (pi() * 3 / 4)))
940  { //LOWER WALL
941  (*w)->SetTransporters1(ch, 0.);
942  (*w)->SetTransporters2(ch, 1.);
943  }
944 
945  else if ((t2 > (-pi() * 3 / 4)) && (t2 <= (-pi() / 4)))
946  {
947  if ((*w)->GetC1() == this)
948  { //LEFT WALL WRT C1
949  (*w)->SetTransporters1(ch, lat);
950  (*w)->SetTransporters2(ch, 0.);
951  }
952  else
953  { //RIGHT WALL WRT C2 (therefore switch transporters 1 and 2)
954  (*w)->SetTransporters1(ch, 0.);
955  (*w)->SetTransporters2(ch, 0.);
956  }
957  }
958  else if ((t2 > (pi() / 4)) && (t2 <= (pi() * 3 / 4))/*t1 > 0.8*/)
959  {
960  if ((*w)->GetC1() == this)
961  { //RIGHT WALL WRT C1
962  (*w)->SetTransporters1(ch, 0.);
963  (*w)->SetTransporters2(ch, 0.);
964  }
965  else
966  { //LEFT WALL WRT C2 (therefore switch transporters 1 and 2)
967  (*w)->SetTransporters1(ch, 0.);
968  (*w)->SetTransporters2(ch, lat);
969  }
970  }
971  }
972  }
973  }
974 }
975 
976 ptree Cell::ToPtree() const
977 {
978  ptree ret = GeometryToPtree();
979  if ( !IsBoundaryPolygon()) {
980  ret.push_back(make_pair("attributes", CellAttributes::ToPtree()));
981  }
982 
983  return ret;
984 }
985 
986 ostream& operator<<(ostream& os, const Cell& c)
987 {
988  c.Print(os);
989  return os;
990 }
991 
992 }
constants.h
Math constants.
std
STL namespace.
SimPT_Sim::Cell
A cell contains walls and nodes.
Definition: Cell.h:48
SimPT_Sim::Node
Node in cell wall.
Definition: Node.h:39
SimPT_Sim::Util
Namespace for miscellaneous utilities.
Definition: PTreeFile.cpp:44
SimPT_Sim::GeoData
Plain data structure with cell geometric data (such as area).
Definition: GeoData.h:35
AreaMoment.h
Interface/Implementation for AreaMoment.
SimPT_Sim::Cell::GetWalls
const std::list< Wall * > & GetWalls() const
Access the cell's walls.
Definition: Cell.h:88
GeoData.h
Interface/Implementation for GeoData.
SimPT_Sim
Namespace for the core simulator.
Definition: CellDivider.cpp:183
SimPT_Sim::Edge
An Edge connects two nodes and is ambidextrous.
Definition: Edge.h:31
Cell.h
Interface for Cell.
SimPT_Sim::Container
Namespace for container related classes.
Definition: CircularIterator.h:27
SimPT_Sim::Cell::GetCircumference
double GetCircumference() const
Return the circumference along the edges.
Definition: Cell.cpp:213
SimPT_Sim::Cell::GeometryToPtree
ptree GeometryToPtree() const
Convert the cell geometry to a ptree.
Definition: Cell.cpp:161
SimPT_Sim::Cell::IsWallNeighbor
bool IsWallNeighbor(Cell *cell) const
Strict neighbor (you're never your own neighbor)
Definition: Cell.cpp:290
SimPT_Sim::Cell::GetSignedArea
double GetSignedArea() const
Return the signed area of the cell.
Definition: Cell.cpp:186
Node.h
Interface for Node.
ConstCircularIterator.h
ConstCircularIterator class and helper functions.
SimPT_Sim::Cell::GetIndex
int GetIndex() const
Return the index.
Definition: Cell.h:76
SimPT_Sim::AreaMoment
Structure with functionality for calculations of the area moment of inertia.
Definition: AreaMoment.h:27
SimPT_Sim::Cell::GetCentroid
std::array< double, 3 > GetCentroid() const
Return the centroid position.
Definition: Cell.cpp:205
Edge.h
Interface for Edge.
SimPT_Sim::CellAttributes::ToPtree
virtual boost::property_tree::ptree ToPtree() const
Convert the cell attributes to a ptree.
Definition: CellAttributes.cpp:103
SimPT_Sim::CellAttributes
Attributes of Cell.
Definition: CellAttributes.h:36
array3.h
Extending std with arithmetic for std::array.
SimPT_Sim::Cell::ToPtree
virtual ptree ToPtree() const
Convert the cell (geometry and attributes) to a ptree.
Definition: Cell.cpp:976
SimPT_Sim::Cell::GetGeoData
GeoData GetGeoData() const
Return GeData (area, centroid, area moment of inertia).
Definition: Cell.cpp:225
SimPT_Sim::Container::next
CircularIterator< T > next(CircularIterator< T > i)
Helper yields the position the iterator would have if moved forward (in circular fashion) by 1 positi...
Definition: CircularIterator.h:94
CellAttributes.h
Interface for CellAttributes.
SimPT_Sim::Cell::GetSumTransporters
double GetSumTransporters(unsigned int ch) const
Sum transporters at this cell's side of the walls.
Definition: Cell.cpp:233
SimPT_Sim::Container::make_const_circular
ConstCircularIterator< typename T::const_iterator > make_const_circular(const T &c)
Helper produces const circular iterator whose range corresponds to the begin and end iterators of a c...
Definition: ConstCircularIterator.h:47
SimPT_Sim::Cell::Move
void Move(const std::array< double, 3 > &a)
Strict neighbor (you're never your own neighbor)
Definition: Cell.cpp:302
SimPT_Sim::Container::prev
CircularIterator< T > prev(CircularIterator< T > i)
Helper yields the position the iterator would have if moved backward (in circular fashion) by 1 posit...
Definition: CircularIterator.h:107
CircularIterator.h
CircularIterator class and helper functions.
SimPT_Sim::Cell::MoveSelfIntersects
bool MoveSelfIntersects(Node *moving_node, std::array< double, 3 > new_pos)
Check for self-intersection when moving_node gets displaced.
Definition: Cell.cpp:309
SimPT_Sim::Cell::GetArea
double GetArea() const
Return the area of the cell.
Definition: Cell.cpp:178
SimPT_Sim::Util::pi
constexpr double pi()
Math constant pi.
Definition: constants.h:29
SimPT_Sim::Wall
A cell wall, runs between cell corner points and consists of wall elements.
Definition: Wall.h:48