AreaFix.cpp
1 #include "AreaFix.hpp"
2 #include <utility>
3 #include <set>
4 
5 namespace
6 {
7  Extensive ComputationalCell2Extensive(ComputationalCell const& cell, double volume, EquationOfState
8  const& eos)
9  {
10  Extensive res;
11  res.mass = cell.density*volume;
12  res.energy = eos.dp2e(cell.density, cell.pressure, cell.tracers)*res.mass +
13  0.5*res.mass*ScalarProd(cell.velocity, cell.velocity);
14  res.momentum = res.mass*cell.velocity;
15  for (size_t i = 0; i < cell.tracers.size(); ++i)
16  res.tracers[i] = res.mass*cell.tracers[i];
17  return res;
18  }
19 
20  int GetEdgeIndex(Tessellation const& tessnew, int n0, int n1, int cell_index)
21  {
22  vector<int> const& edges = tessnew.GetCellEdges(cell_index);
23  for (size_t j = 0; j<edges.size(); ++j)
24  {
25  Edge const& etemp = tessnew.GetEdge(edges[j]);
26  if ((tessnew.GetOriginalIndex(etemp.neighbors.first) == n0&&
27  tessnew.GetOriginalIndex(etemp.neighbors.second) == n1) ||
28  (tessnew.GetOriginalIndex(etemp.neighbors.first) == n1&&
29  tessnew.GetOriginalIndex(etemp.neighbors.second) == n0))
30  return edges[j];
31  }
32  return -1;
33  }
34 
35  Vector2D FixPeriodicLeap(Vector2D const& point, Vector2D const& vel, double dt, OuterBoundary const& outer)
36  {
37  Vector2D res(0, 0);
38  if ((point.x + vel.x*dt)>outer.GetGridBoundary(Right))
39  res.x = outer.GetGridBoundary(Right) - outer.GetGridBoundary(Left);
40  if ((point.x + vel.x*dt) < outer.GetGridBoundary(Left))
41  res.x = outer.GetGridBoundary(Left) - outer.GetGridBoundary(Right);
42  if ((point.y + vel.y*dt) > outer.GetGridBoundary(Up))
43  res.y = outer.GetGridBoundary(Up) - outer.GetGridBoundary(Down);
44  if ((point.y + vel.y*dt) < outer.GetGridBoundary(Down))
45  res.y = outer.GetGridBoundary(Down) - outer.GetGridBoundary(Up);
46  return res;
47  }
48 
49  void GetAdjacentNeigh(Tessellation const& tessold, Edge const& edge, int cell_index,
50  Vector2D const& added, int &n0other, int &n1other)
51  {
52  vector<int> const& edges = tessold.GetCellEdges(cell_index);
53  vector<double> dist;
54  Vector2D point;
55  dist.reserve(edges.size());
56  for (size_t i = 0; i < edges.size(); ++i)
57  {
58  Edge const& e = tessold.GetEdge(edges[i]);
59  if (e.neighbors.first == edge.neighbors.first&&e.neighbors.second == edge.neighbors.second)
60  dist.push_back(edge.GetLength() * 1000);
61  else
62  dist.push_back(min(min(min(abs(e.vertices.first - edge.vertices.first + added),
63  abs(e.vertices.second - edge.vertices.first + added)),
64  abs(e.vertices.first - edge.vertices.second + added)),
65  abs(e.vertices.second - edge.vertices.second + added)));
66  }
67  vector<size_t> indeces = sort_index(dist);
68  // The location of the edges to the left and the right of the old edge
69  n0other = (tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[0]]).neighbors.first) == cell_index) ?
70  tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[0]]).neighbors.second) :
71  tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[0]]).neighbors.first);
72  n1other = (tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[1]]).neighbors.first) == cell_index) ?
73  tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[1]]).neighbors.second) :
74  tessold.GetOriginalIndex(tessold.GetEdge(edges[indeces[1]]).neighbors.first);
75  }
76 
77  int NewEdgeIndex(Tessellation const& tessold, Tessellation const& tessnew, int cell_index,
78  Edge const& edge, int other_index)
79  {
80  int n0, n1;
81  int res = -1;
82  bool rigid = tessold.GetOriginalIndex(edge.neighbors.first) == tessold.GetOriginalIndex(edge.neighbors.second) ?
83  true : false;
84  // Find n0 and n1 the new flip neighbors
85  GetAdjacentNeigh(tessold, edge, cell_index, Vector2D(0, 0), n0, n1);
86  // Do we have two neighboring rigid edges?
87  if (rigid && (n0 == cell_index || n1 == cell_index))
88  return -1;
89 
90  // Are we flipping a real edge to a rigid one?
91  if (n0 == cell_index)
92  {
93  n0 = n1;
94  n1 = cell_index;
95  }
96  if (n1 == cell_index)
97  {
98  vector<int> const& newedges = tessnew.GetCellEdges(n0);
99  for (size_t i = 0; i < newedges.size(); ++i)
100  {
101  Edge const& newedge = tessnew.GetEdge(newedges[i]);
102  if (tessnew.GetOriginalIndex(newedge.neighbors.first) ==
103  tessnew.GetOriginalIndex(newedge.neighbors.second))
104  return newedges[i];
105  }
106  return -1;
107  }
108 
109  vector<int> const& edges = tessnew.GetCellEdges(n0);
110  for (size_t i = 0; i < edges.size(); ++i)
111  {
112  Edge const& edge_temp = tessnew.GetEdge(edges[i]);
113  int n0new = tessnew.GetOriginalIndex(edge_temp.neighbors.first), n1new = tessnew.GetOriginalIndex(edge_temp.neighbors.second);
114  if (n0new == n1 || n1new == n1)
115  res = edges[i];
116  }
117  if (res == -1)
118  return -1;
119  if (rigid)
120  return res;
121  int n0new, n1new;
122  GetAdjacentNeigh(tessnew, tessnew.GetEdge(res), n0, Vector2D(0, 0), n0new, n1new);
123  if ((n0new == cell_index&&n1new == other_index) || (n1new == cell_index&&other_index == n0new))
124  return res;
125  // We didn't find a new edge
126  return -1;
127  }
128 
129  bool FirstVertice(Edge edgenew, Tessellation const& tessnew, int cell, Vector2D const& edge_added)
130  {
131  edgenew.vertices.first += edge_added;
132  edgenew.vertices.second += edge_added;
133  const double R = 1e-8*tessnew.GetWidth(cell);
134  vector<int> const& edges = tessnew.GetCellEdges(cell);
135  for (size_t i = 0; i < edges.size(); ++i)
136  {
137  Edge edge = tessnew.GetEdge(edges[i]);
138  if (edge.vertices.first.distance(edgenew.vertices.first) < R)
139  return true;
140  if (edge.vertices.first.distance(edgenew.vertices.second) < R)
141  return false;
142  if (edge.vertices.second.distance(edgenew.vertices.first) < R)
143  return true;
144  if (edge.vertices.second.distance(edgenew.vertices.second) < R)
145  return false;
146  }
147  throw UniversalError("Couldn't find first vertice");
148  }
149 
150  Edge FindOtherEdge(Tessellation const& tess, int cell_index, int other)
151  {
152  vector<int> const& edges = tess.GetCellEdges(cell_index);
153  for (size_t i = 0; i < edges.size(); ++i)
154  {
155  Edge const& edge = tess.GetEdge(edges[i]);
156  if (tess.GetOriginalIndex(edge.neighbors.first) == other)
157  return edge;
158  if (tess.GetOriginalIndex(edge.neighbors.second) == other)
159  return edge;
160  }
161  throw UniversalError("Couldn't find other edge");
162  }
163 
164  std::pair<Vector2D, Vector2D> TrianglesArea(Edge const& eold, Edge const &enew, Tessellation const& tessnew,
165  Tessellation const& tessold, Vector2D const& new_edge_added, int cell_index, Vector2D const& cell_index_added)
166  { // first is the change in the old and the second the change in the new
167  // x is e1.first y is e1.second
168  std::pair<Vector2D, Vector2D> res;
169  boost::array<Vector2D, 4> points;
170  double area_scale = 0;
171  double sum = 0;
172  bool rigid = tessold.GetOriginalIndex(eold.neighbors.first) ==
173  tessold.GetOriginalIndex(eold.neighbors.second);
174  Vector2D oldmeshpoint = (tessold.GetPointNo() < eold.neighbors.first) && rigid ?
175  tessold.GetMeshPoint(eold.neighbors.second) : tessold.GetMeshPoint(eold.neighbors.first);
176  TripleConstRef<Vector2D> temp(eold.vertices.first, eold.vertices.second, oldmeshpoint);
177  if (orient2d(temp) > 0)
178  {
179  points[0] = eold.vertices.first;
180  points[1] = eold.vertices.second;
181  }
182  else
183  {
184  points[1] = eold.vertices.first;
185  points[0] = eold.vertices.second;
186  }
187  Vector2D newmeshpoint = tessnew.GetMeshPoint(enew.neighbors.first);
188  TripleConstRef<Vector2D> temp2(enew.vertices.first, enew.vertices.second, newmeshpoint);
189  if (orient2d(temp2) > 0)
190  {
191  points[3] = enew.vertices.first + new_edge_added;
192  points[2] = enew.vertices.second + new_edge_added;
193  }
194  else
195  {
196  points[2] = enew.vertices.first + new_edge_added;
197  points[3] = enew.vertices.second + new_edge_added;
198  }
199  int npoints = tessold.GetPointNo();
200  Vector2D toadd = new_edge_added - cell_index_added;
201  bool first = FirstVertice(enew, tessnew, cell_index, toadd);
202  if (eold.neighbors.first > npoints)
203  {
204  first = !first;
205  if (rigid)
206  {
207  Vector2D tempp = points[0];
208  points[0] = points[1];
209  points[1] = tempp;
210  }
211  }
212  Vector2D third_point = (first ? enew.vertices.first : enew.vertices.second) + new_edge_added;
213  {
214  TripleConstRef<Vector2D> temp3 = TripleConstRef<Vector2D>(points[0], third_point, points[1]);
215  res.first.x = 0.5*orient2d(temp3);
216  area_scale = std::abs(res.first.x);
217  sum = res.first.x;
218  }
219  {
220  third_point = (first ? enew.vertices.second : enew.vertices.first) + new_edge_added;
221  TripleConstRef<Vector2D> temp4 = TripleConstRef<Vector2D>(points[1], third_point, points[0]);
222  res.first.y = 0.5*orient2d(temp4);
223  area_scale = std::max(area_scale, std::abs(res.first.y));
224  sum += res.first.y;
225  }
226  bool newrigid = tessnew.GetOriginalIndex(enew.neighbors.first) == tessnew.GetOriginalIndex(enew.neighbors.second);
227  int other = tessnew.GetOriginalIndex(enew.neighbors.first);
228  Edge otheredge = FindOtherEdge(tessold, cell_index, other);
229  Vector2D added;
230  if (tessold.GetOriginalIndex(otheredge.neighbors.first) == other)
231  added = tessold.GetMeshPoint(other) - tessold.GetMeshPoint(otheredge.neighbors.first);
232  else
233  added = tessold.GetMeshPoint(other) - tessold.GetMeshPoint(otheredge.neighbors.second);
234  first = FirstVertice(eold, tessold, other, added);
235 
236  if (newrigid&&enew.neighbors.first > npoints)
237  first = !first;
238  if (enew.neighbors.first < npoints || (tessnew.GetOriginalIndex(enew.neighbors.first) !=
239  tessnew.GetOriginalIndex(enew.neighbors.second)))
240  {
241  third_point = first ? eold.vertices.first : eold.vertices.second;
242  TripleConstRef<Vector2D> temp5(points[2], third_point, points[3]);
243  res.second.x = 0.5*orient2d(temp5);
244  area_scale = std::max(area_scale, std::abs(res.second.x));
245  sum += res.second.x;
246  }
247 
248  if (enew.neighbors.second < npoints || (tessnew.GetOriginalIndex(enew.neighbors.first) !=
249  tessnew.GetOriginalIndex(enew.neighbors.second)))
250  {
251  third_point = !first ? eold.vertices.first : eold.vertices.second;
252  TripleConstRef<Vector2D> temp6(points[3], third_point, points[2]);
253  res.second.y = 0.5*orient2d(temp6);
254  area_scale = std::max(area_scale, std::abs(res.second.y));
255  sum += res.second.y;
256  }
257  if (std::abs(sum) > 1e-5*area_scale)
258  {
259  UniversalError eo("Bad total area in TrianglesArea");
260  eo.AddEntry("area_scale", area_scale);
261  eo.AddEntry("sum", sum);
262  eo.AddEntry("old neigh 0", eold.neighbors.first);
263  eo.AddEntry("old neigh 1", eold.neighbors.second);
264  throw eo;
265  }
266  return res;
267  }
268 
269  Vector2D FixNewEdgeLeap(Tessellation const& tessold, int real_cell, int cell_index)
270  {
271  vector<int> const& edges = tessold.GetCellEdges(cell_index);
272  for (size_t i = 0; i < edges.size(); ++i)
273  {
274  Edge const& edge = tessold.GetEdge(edges[i]);
275  if (tessold.GetOriginalIndex(edge.neighbors.first) == real_cell)
276  return tessold.GetMeshPoint(edge.neighbors.first) - tessold.GetMeshPoint(real_cell);
277  if (tessold.GetOriginalIndex(edge.neighbors.second) == real_cell)
278  return tessold.GetMeshPoint(edge.neighbors.second) - tessold.GetMeshPoint(real_cell);
279  }
280  throw UniversalError("Couldn't find New edge leap");
281  }
282 
283  void GetCellIndex(Edge const& edge, Tessellation const& tess, int &cell_index, int &other_index)
284  {
285  int npoints = tess.GetPointNo();
286  if (edge.neighbors.first < npoints)
287  {
288  cell_index = edge.neighbors.first;
289  other_index = tess.GetOriginalIndex(edge.neighbors.second);
290  }
291  else
292  {
293  cell_index = edge.neighbors.second;
294  other_index = tess.GetOriginalIndex(edge.neighbors.first);
295  }
296  }
297 
298  double GetdAFlux(int mid_index, Tessellation const& tessmid, double dt,
299  int other_index, vector<Vector2D> const& fv)
300  {
301  double res = 0;
302  if (mid_index >= 0)
303  {
304  Edge const& edge2 = tessmid.GetEdge(mid_index);
305  const Vector2D norm = tessmid.GetMeshPoint(edge2.neighbors.second) - tessmid.GetMeshPoint(edge2.neighbors.first);
306  res = ScalarProd(norm, fv[static_cast<size_t>(mid_index)])*dt*edge2.GetLength() / abs(norm);
307  if (tessmid.GetOriginalIndex(edge2.neighbors.first) == other_index)
308  res *= -1;
309  }
310  return res;
311  }
312 
313  void AreaFixEdgeDisappear(Tessellation const& tessold, Tessellation const& tessnew, int cell_index, int other_index,
314  Edge const& edge, OuterBoundary const& outer, int npoints, vector<Vector2D> const& pointvelocity, double dt,
315  double dA_flux, int mid_index, Tessellation const& tessmid, vector<Vector2D> const& fv,
316  vector<ComputationalCell> const& cells, EquationOfState const& eos, vector<Extensive> &res,
317  std::set<std::pair<int, int > > &flipped_set, Vector2D const& cell_added)
318  {
319  // Is there a corresponding new edge? An edge flip
320  int new_edge = NewEdgeIndex(tessold, tessnew, cell_index, edge, other_index);
321  if (new_edge < 0)
322  return;
323  bool both_real = edge.neighbors.first < npoints && edge.neighbors.second < npoints;
324  std::pair<int, int> neighbors(tessold.GetOriginalIndex(edge.neighbors.first), tessold.GetOriginalIndex(edge.neighbors.second));
325  if (flipped_set.count(neighbors) > 0)
326  return;
327  Edge edge_new = tessnew.GetEdge(new_edge);
328  // Did it jump?
329  Vector2D addedother;
330  if (outer.GetBoundaryType() == Periodic)
331  {
332  if (edge_new.neighbors.first < npoints)
333  {
334  addedother = -1 * FixPeriodicLeap(tessnew.GetMeshPoint(edge_new.neighbors.first),
335  pointvelocity[static_cast<size_t>(edge_new.neighbors.first)], -dt, outer);
336  addedother += FixNewEdgeLeap(tessold, edge_new.neighbors.first, cell_index);
337  }
338  else
339  {
340  addedother = -1 * FixPeriodicLeap(tessnew.GetMeshPoint(edge_new.neighbors.second),
341  pointvelocity[static_cast<size_t>(edge_new.neighbors.second)], -dt, outer);
342  addedother += FixNewEdgeLeap(tessold, edge_new.neighbors.second, cell_index);
343  }
344  }
345  std::pair<Vector2D, Vector2D> areas = TrianglesArea(edge, edge_new, tessnew, tessold, addedother, cell_index,
346  cell_added);
347  vector<double> dA(4);
348  if (mid_index >= 0 && tessmid.GetOriginalIndex(tessmid.GetEdge(mid_index).neighbors.first) == other_index)
349  dA_flux *= -1;
350  if (mid_index >= 0 && tessmid.GetOriginalIndex(tessmid.GetEdge(mid_index).neighbors.first) == tessold.GetOriginalIndex(edge.neighbors.second))
351  dA_flux *= -1;
352  dA[0] = areas.first.x - dA_flux;
353  dA[1] = areas.first.y + dA_flux;
354  dA[2] = areas.second.x;
355  dA[3] = areas.second.y;
356 
357  // Did the new edge exist in the mid time step? This should mean mid_index=-1
358  int mid_index2 = GetEdgeIndex(tessmid, tessnew.GetOriginalIndex(edge_new.neighbors.first),
359  tessnew.GetOriginalIndex(edge_new.neighbors.second), tessnew.GetOriginalIndex(edge_new.neighbors.first));
360 
361  if (tessnew.GetOriginalIndex(edge_new.neighbors.first) == tessnew.GetOriginalIndex(edge_new.neighbors.second))
362  mid_index2 = -1;
363 
364  double dA_flux2 = GetdAFlux(mid_index2, tessmid, dt, tessnew.GetOriginalIndex(edge_new.neighbors.first), fv);
365  dA[2] += dA_flux2;
366  dA[3] -= dA_flux2;
367 
368  // The conserved to remove
369  Extensive TotalRemoved(cells[0].tracers);
370  double total_added_area = 0;
371  if (dA[0] < 0)
372  {
373  Extensive toremove = ComputationalCell2Extensive(cells[
374  static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.first))], -dA[0], eos);
375  res[static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.first))] -= toremove;
376  TotalRemoved += toremove;
377  }
378  else
379  total_added_area += dA[0];
380  if (dA[1] < 0)
381  {
382  Extensive toremove = ComputationalCell2Extensive(cells[
383  static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.second))], -dA[1], eos);
384  res[static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.second))] -= toremove;
385  TotalRemoved += toremove;
386  }
387  else
388  total_added_area += dA[1];
389 
390  if (dA[2] < 0)
391  {
392  Extensive toremove = ComputationalCell2Extensive(cells[
393  static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.first))], -dA[2], eos);
394  res[static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.first))] -= toremove;
395  TotalRemoved += toremove;
396  }
397  else
398  total_added_area += dA[2];
399  if (dA[3] < 0)
400  {
401  Extensive toremove = ComputationalCell2Extensive(cells[
402  static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.second))], -dA[3], eos);
403  res[static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.second))] -= toremove;
404  TotalRemoved += toremove;
405  }
406  else
407  total_added_area += dA[3];
408 
409  // Add the conserved
410  if (dA[0] > 0)
411  res[static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.first))] += (dA[0] / total_added_area)*TotalRemoved;
412  if (dA[1] > 0)
413  res[static_cast<size_t>(tessold.GetOriginalIndex(edge.neighbors.second))] += (dA[1] / total_added_area)*TotalRemoved;
414  if (dA[2] > 0)
415  if ((edge_new.neighbors.first < tessnew.GetPointNo()) ||
416  (tessnew.GetOriginalIndex(edge_new.neighbors.first) != tessnew.GetOriginalIndex(edge_new.neighbors.second)))
417  res[static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.first))] += (dA[2] / total_added_area)*TotalRemoved;
418  if (dA[3] > 0)
419  if ((edge_new.neighbors.second < tessnew.GetPointNo()) ||
420  (tessnew.GetOriginalIndex(edge_new.neighbors.first) != tessnew.GetOriginalIndex(edge_new.neighbors.second)))
421  res[static_cast<size_t>(tessnew.GetOriginalIndex(edge_new.neighbors.second))] += (dA[3] / total_added_area)*TotalRemoved;
422  if (!both_real)
423  {
424  flipped_set.insert(neighbors);
425  std::pair<int, int> neigh_temp(neighbors.second, neighbors.first);
426  flipped_set.insert(neigh_temp);
427  }
428  }
429 
430  double EdgesArea(Edge const& e1, Edge const &e2, Vector2D const& pointold, Vector2D const& pointnew)
431  {
432  // e1 is old tess, e2 is new tess
433  boost::array<Vector2D, 4> edge_points;
434  double res = 0;
435  TripleConstRef<Vector2D> temp1(pointold, e1.vertices.first, e1.vertices.second);
436  TripleConstRef<Vector2D> temp2(pointnew, e2.vertices.first, e2.vertices.second);
437  if (orient2d(temp1) > 0)
438  {
439  edge_points[0] = e1.vertices.first;
440  edge_points[1] = e1.vertices.second;
441  }
442  else
443  {
444  edge_points[1] = e1.vertices.first;
445  edge_points[0] = e1.vertices.second;
446  }
447  if (orient2d(temp2) > 0)
448  {
449  edge_points[2] = e2.vertices.first;
450  edge_points[3] = e2.vertices.second;
451  }
452  else
453  {
454  edge_points[3] = e2.vertices.first;
455  edge_points[2] = e2.vertices.second;
456  }
457  TripleConstRef<Vector2D> tri1(edge_points[0], edge_points[3], edge_points[1]);
458  TripleConstRef<Vector2D> tri2(edge_points[0], edge_points[2], edge_points[3]);
459  res += orient2d(tri1);
460  res += orient2d(tri2);
461  return res*0.5;
462  }
463 
464  ComputationalCell GetDonorCell(ComputationalCell const& c0, ComputationalCell const& c1,
465  bool seconddonor, EquationOfState const& /*eos*/)
466  {
467  ComputationalCell p_mid;
468  double ratio = c1.density / c0.density;
469  if (ratio<2 && ratio>0.5)
470  p_mid.density = 0.5*(c1.density + c0.density);
471  else
472  {
473  if (seconddonor)
474  p_mid.density = c1.density;
475  else
476  p_mid.density = c0.density;
477  }
478  ratio = c1.pressure / c0.pressure;
479  if (ratio<2 && ratio>0.5)
480  p_mid.pressure = 0.5*(c1.pressure + c0.pressure);
481  else
482  {
483  if (seconddonor)
484  p_mid.pressure = c1.pressure;
485  else
486  p_mid.pressure = c0.pressure;
487  }
488  ratio = c1.velocity.x / c0.velocity.x;
489  if (ratio<2 && ratio>0.5)
490  p_mid.velocity.x = 0.5*(c1.velocity.x + c0.velocity.x);
491  else
492  {
493  if (seconddonor)
494  p_mid.velocity.x = c1.velocity.x;
495  else
496  p_mid.velocity.x = c0.velocity.x;
497  }
498  ratio = c1.velocity.y / c0.velocity.y;
499  if (ratio<2 && ratio>0.5)
500  p_mid.velocity.y = 0.5*(c1.velocity.y + c0.velocity.y);
501  else
502  {
503  if (seconddonor)
504  p_mid.velocity.y = c1.velocity.y;
505  else
506  p_mid.velocity.y = c0.velocity.y;
507  }
508  for (size_t i = 0; i < c0.tracers.size(); ++i)
509  {
510  ratio = c0.tracers[i] / c1.tracers[i];
511  if (ratio<2 && ratio>0.5)
512  p_mid.tracers[i] = 0.5*c0.tracers[i] * (1 + 1 / ratio);
513  else
514  if (seconddonor)
515  p_mid.tracers[i] = c1.tracers[i];
516  else
517  p_mid.tracers[i] = c0.tracers[i];
518  }
519  return p_mid;
520  }
521 }
522 
523 vector<Extensive> FluxFix2(Tessellation const& tessold, Tessellation const& tessmid,
524  Tessellation const& tessnew, vector<Vector2D> const& pointvelocity, double dt,
525  vector<ComputationalCell> const& cells, vector<Extensive> const& /*fluxes*/,
526  vector<Vector2D> const& fv, OuterBoundary const& outer, EquationOfState const& eos)
527 {
528  size_t npoints = static_cast<size_t>(tessold.GetPointNo());
529  vector<Extensive> res(static_cast<size_t>(npoints), Extensive(cells[0].tracers));
530  // Fix the fluxes
531  Vector2D temp0, temp1;
532  std::set<std::pair<int, int > > flipped_set, only_mid;
533  size_t nedgesold = static_cast<size_t>(tessold.GetTotalSidesNumber());
534  for (size_t i = 0; i < nedgesold; ++i)
535  {
536  Edge const& edge = tessold.GetEdge(static_cast<int>(i));
537 
538  int cell_index, other_index;
539  GetCellIndex(edge, tessold, cell_index, other_index);
540  int mid_index = GetEdgeIndex(tessmid, cell_index, other_index, cell_index);
541  int new_index = GetEdgeIndex(tessnew, cell_index, other_index, cell_index);
542  double dA_flux = GetdAFlux(mid_index, tessmid, dt, other_index, fv);
543 
544  Vector2D real_p;
545  real_p = tessold.GetMeshPoint(cell_index);
546  Vector2D added(0, 0);
547  // check if periodic jumped
548  if (outer.GetBoundaryType() == Periodic)
549  added = FixPeriodicLeap(real_p, pointvelocity[static_cast<size_t>(cell_index)], dt, outer);
550 
551  // Did the edge disappear? An edge flip
552  if (new_index < 0)
553  {
554  AreaFixEdgeDisappear(tessold, tessnew, cell_index, other_index, edge, outer, static_cast<int>(npoints), pointvelocity, dt,
555  dA_flux, mid_index, tessmid, fv, cells, eos, res, flipped_set, added);
556  continue;
557  }
558  Vector2D real_p_new = tessnew.GetMeshPoint(cell_index);
559  Edge edge_other = tessnew.GetEdge(new_index);
560  if (outer.GetBoundaryType() == Periodic)
561  {
562  edge_other.vertices.first += added;
563  edge_other.vertices.second += added;
564  real_p_new += added;
565  }
566  double area = EdgesArea(edge, edge_other, real_p, real_p_new);
567  ComputationalCell p_mid = GetDonorCell(cells[static_cast<size_t>(cell_index)], cells[static_cast<size_t>(other_index)],
568  (area - dA_flux)>0, eos);
569  Extensive toadd = ComputationalCell2Extensive(p_mid, (-dA_flux + area), eos);
570 
571  if (other_index == edge.neighbors.first || other_index == edge.neighbors.second)
572  res[static_cast<size_t>(other_index)] -= toadd;
573  if (cell_index == edge.neighbors.first || cell_index == edge.neighbors.second)
574  res[static_cast<size_t>(cell_index)] += toadd;
575  if (mid_index < 0)
576  {
577  // Was there a flip in mid step?
578  int mid_edge = NewEdgeIndex(tessold, tessmid, cell_index, edge, other_index);
579  if (mid_edge < 0)
580  continue;
581  Edge const& edge2 = tessmid.GetEdge(mid_edge);
582  if (tessmid.GetOriginalIndex(edge2.neighbors.first) == tessmid.GetOriginalIndex(edge2.neighbors.second))
583  continue;
584  double dA_flux2 = GetdAFlux(mid_edge, tessmid, dt, tessmid.GetOriginalIndex(edge2.neighbors.second), fv);
585  Extensive toadd2(cells[0].tracers);
586  if (dA_flux2 > 0)
587  toadd2 = ComputationalCell2Extensive(cells[static_cast<size_t>(tessmid.GetOriginalIndex(edge2.neighbors.second))],
588  dA_flux2, eos);
589  else
590  toadd2 = ComputationalCell2Extensive(cells[static_cast<size_t>(tessmid.GetOriginalIndex(edge2.neighbors.first))],
591  dA_flux2, eos);
592  std::pair<int, int> mid_neigh(tessmid.GetOriginalIndex(edge2.neighbors.first), tessmid.GetOriginalIndex(edge2.neighbors.second));
593  std::pair<int, int> mid_neigh2(mid_neigh.second, mid_neigh.first);
594  if (only_mid.count(mid_neigh) == 0)
595  {
596  only_mid.insert(mid_neigh);
597  only_mid.insert(mid_neigh2);
598  res[static_cast<size_t>(tessmid.GetOriginalIndex(edge2.neighbors.second))] += toadd2;
599  res[static_cast<size_t>(tessmid.GetOriginalIndex(edge2.neighbors.first))] -= toadd2;
600  }
601  }
602  }
603  return res;
604 }
Extensive
Extensive variables.
Definition: extensive.hpp:18
sort_index
void sort_index(const vector< T > &arr, vector< int > &res)
Returns the indeces of a sort.
Definition: utils.hpp:500
Tessellation
Abstract class for tessellation.
Definition: tessellation.hpp:21
Tessellation::GetPointNo
virtual int GetPointNo(void) const =0
Get Total number of mesh generating points.
AreaFix.hpp
Fixes the area inconsistency problem.
Tessellation::GetOriginalIndex
virtual int GetOriginalIndex(int point) const
Returns the original index of the duplicated point.
Definition: tessellation.cpp:5
Extensive::momentum
Vector2D momentum
momentum, in relativity it is = rho*h*gamma*v
Definition: extensive.hpp:31
UniversalError
Container for error reports.
Definition: universal_error.hpp:14
Extensive::mass
double mass
rest mass times gamma
Definition: extensive.hpp:25
Edge
Interface between two cells.
Definition: Edge.hpp:13
Tessellation::GetMeshPoint
virtual Vector2D GetMeshPoint(int index) const =0
Returns Position of mesh generating point.
max
double max(vector< double > const &v)
returns the maximal term in a vector
Definition: utils.cpp:52
Extensive::energy
double energy
energy, in relativity it is = rho*h*gamma^2-P-rho
Definition: extensive.hpp:28
Vector2D::y
double y
Component in the y direction.
Definition: geometry.hpp:48
TripleConstRef
A collection of three identical references.
Definition: triplet.hpp:12
ComputationalCell::tracers
tvector tracers
Tracers (can transfer from one cell to another)
Definition: computational_cell_2d.hpp:40
Tessellation::GetEdge
virtual Edge const & GetEdge(int index) const =0
Returns edge (interface between cells)
ComputationalCell::pressure
double pressure
Pressure.
Definition: computational_cell_2d.hpp:34
Extensive::tracers
tvector tracers
tracers
Definition: extensive.hpp:34
ScalarProd
double ScalarProd(Vector3D const &v1, Vector3D const &v2)
Scalar product of two vectors.
Definition: Vector3D.cpp:185
OuterBoundary::GetBoundaryType
virtual BoundaryType GetBoundaryType(void) const =0
Returns the boundary type.
EquationOfState
Base class for equation of state.
Definition: equation_of_state.hpp:20
Edge::vertices
std::pair< Vector2D, Vector2D > vertices
Points at the ends of the edge.
Definition: Edge.hpp:18
Tessellation::GetWidth
virtual double GetWidth(int index) const =0
Returns the effective width of a cell.
min
double min(vector< double > const &v)
Returns the minimal term in a vector.
Definition: utils.cpp:44
Tessellation::GetTotalSidesNumber
virtual int GetTotalSidesNumber(void) const =0
Returns the total number of faces.
Edge::neighbors
std::pair< int, int > neighbors
Neighboring cells.
Definition: Edge.hpp:21
EquationOfState::dp2e
virtual double dp2e(double d, double p, tvector const &tracers=tvector(), vector< string > const &tracernames=vector< string >()) const =0
Calculates the thermal energy per unit mass.
TripleConstRef::first
const T & first
Reference to first item.
Definition: triplet.hpp:17
abs
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
Tessellation::GetCellEdges
virtual vector< int > const & GetCellEdges(int index) const =0
Returns the indexes of a cell&#39;s edges.
ComputationalCell::velocity
Vector2D velocity
Velocity.
Definition: computational_cell_2d.hpp:37
OuterBoundary
Abstract class for geometric boundary conditions for the tessellation.
Definition: OuterBoundary.hpp:19
Vector2D
2D Mathematical vector
Definition: geometry.hpp:15
Edge::GetLength
double GetLength(void) const
Returns the length of the edge.
Definition: Edge.cpp:26
OuterBoundary::GetGridBoundary
virtual double GetGridBoundary(Directions dir) const =0
Returns the boundary coordinate.
orient2d
double orient2d(const TripleConstRef< Vector2D > &points)
Checks whether 3 given points are on a counterclockwise circle, clockwise circle or colinear...
Definition: geotests.cpp:76
FluxFix2
vector< Extensive > FluxFix2(Tessellation const &tessold, Tessellation const &tessmid, Tessellation const &tessnew, vector< Vector2D > const &pointvelocity, double dt, vector< ComputationalCell > const &cells, vector< Extensive > const &fluxes, vector< Vector2D > const &fv, OuterBoundary const &outer, EquationOfState const &eos)
Fixes the flux to propely converge second order.
Definition: AreaFix.cpp:523
Vector2D::x
double x
Component in the x direction.
Definition: geometry.hpp:45
ComputationalCell
Computational cell.
Definition: computational_cell_2d.hpp:22
Generated on Thu Jun 4 2020 18:12:16 by   doxygen 1.8.13