openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
SF_fem_utils_emi.h
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
29 #if WITH_EMI_MODEL
30 
31 #ifndef _SF_FEM_EMI_H
32 #define _SF_FEM_EMI_H
33 
34 #include "SF_container.h"
35 #include "SF_fem_utils.h"
36 #include "SF_init.h"
37 namespace SF {
38 
49 inline double computeTriangleArea(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3)
50 {
51  // Define two edges of the triangle
52  Point edge1 = p2 - p1;
53  Point edge2 = p3 - p1;
54 
55  // Compute the cross product
56  Point crossProduct = cross(edge1,edge2);
57 
58  // Compute the magnitude of the cross product
59  double crossProductMagnitude = mag(crossProduct);
60 
61  return 0.5 * crossProductMagnitude;
62 }
63 
74 inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& P, vector<double>& interpolationCoefficients, double & lengthSquared)
75 {
76  SF::Point d = p2 - p1;
77  lengthSquared = inner_prod(d, d);
78 
79  SF::Point v = P - p1;
80  double lambda2 = inner_prod(v, d) / lengthSquared;
81  double lambda1 = 1.0 - lambda2;
82 
83  interpolationCoefficients[0] = lambda1;
84  interpolationCoefficients[1] = lambda2;
85 }
86 
98 inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3, const SF::Point& P, vector<double>& interpolationCoefficients, double & area)
99 {
100 
101  area = computeTriangleArea(p1, p2, p3);
102 
103  // Vectors for the edges of the triangle
104  SF::Point v0 = p2 - p1;
105  SF::Point v1 = p3 - p1;
106  SF::Point v2 = P - p1;
107 
108  // Dot products
109  double d00 = inner_prod(v0,v0);
110  double d01 = inner_prod(v0,v1);
111  double d11 = inner_prod(v1,v1);
112  double d20 = inner_prod(v2,v0);
113  double d21 = inner_prod(v2,v1);
114 
115  // Compute barycentric coordinates
116  double denom = d00 * d11 - d01 * d01;
117  double beta = (d11 * d20 - d01 * d21) / denom;
118  double gamma = (d00 * d21 - d01 * d20) / denom;
119  double alpha = 1.0 - beta - gamma;
120  interpolationCoefficients[0] = alpha; // α
121  interpolationCoefficients[1] = beta; // β
122  interpolationCoefficients[2] = gamma; // ζ
123 }
124 
125 
138 inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3, const SF::Point& p4, const SF::Point& P, vector<double>& interpolationCoefficients, double & area)
139 {
140  // Area: compute as sum of two triangles
141  double area1 = computeTriangleArea(p1, p2, p4);
142  double area2 = computeTriangleArea(p2, p3, p4);
143  area = area1 + area2;
144 
145  // Define local axes
146  SF::Point origin = p1;
147  SF::Point u_dir = p2 - p1;
148  SF::Point v_dir = p4 - p1;
149 
150  // Project P into the local (u,v) frame
151  SF::Point d = P - origin;
152  double u = inner_prod(d, u_dir) / inner_prod(u_dir, u_dir);
153  double v = inner_prod(d, v_dir) / inner_prod(v_dir, v_dir);
154 
155  // Clamp u and v to [0,1] to avoid extrapolation (optional)
156  u = SF::clamp(u, 0.0, 1.0);
157  v = SF::clamp(v, 0.0, 1.0);
158 
159  // Bilinear shape functions
160  interpolationCoefficients[0] = (1 - u) * (1 - v); // α
161  interpolationCoefficients[1] = u * (1 - v); // β
162  interpolationCoefficients[2] = u * v; // γ
163  interpolationCoefficients[3] = (1 - u) * v; // ζ
164 }
165 
176 template<class S>
177 inline void compute_integrate_matrix_barycentric(vector<SF::Point> face_coordinates, SF_int nnodes, dmat<SF_real> & ebuff, dmat<SF_real> & ebuff_s, dmat<SF_real> & ebuff_counter, S mass_scale)
178 {
179  //compute the barycentric coordinates: calculateCentroid
180  double x = 0;
181  double y = 0;
182  double z = 0;
183  for (int i = 0; i < nnodes; ++i)
184  {
185  x += face_coordinates[i].x;
186  y += face_coordinates[i].y;
187  z += face_coordinates[i].z;
188  }
189  SF::Point b; b.x = x/nnodes; b.y = y/nnodes; b.z = z/nnodes;
190 
191  vector<double> interpolationCoefficients(nnodes);
192  double area, weighted_area;
193  if(nnodes==2){
194  compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], b, interpolationCoefficients, area);
195  }
196  else if(nnodes==3)
197  {
198  compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], face_coordinates[2], b, interpolationCoefficients, area);
199  }
200  else if(nnodes==4){
201  compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], face_coordinates[2], face_coordinates[3], b, interpolationCoefficients, area);
202  }
203 
204  weighted_area = area/nnodes;
205 
206  if(nnodes==2){
207  ebuff[0][0] = interpolationCoefficients[0];
208  ebuff[0][1] = interpolationCoefficients[1];
209 
210  ebuff_counter[0][0] = interpolationCoefficients[0];
211  ebuff_counter[0][1] = interpolationCoefficients[1];
212 
213  ebuff_s[0][0] = mass_scale*weighted_area;
214  ebuff_s[0][1] = mass_scale*weighted_area;
215  }
216  else if(nnodes==3){
217  ebuff[0][0] = interpolationCoefficients[0];
218  ebuff[0][1] = interpolationCoefficients[1];
219  ebuff[0][2] = interpolationCoefficients[2];
220 
221  ebuff_counter[0][0] = interpolationCoefficients[0];
222  ebuff_counter[0][1] = interpolationCoefficients[1];
223  ebuff_counter[0][2] = interpolationCoefficients[2];
224 
225  ebuff_s[0][0] = mass_scale*weighted_area;
226  ebuff_s[0][1] = mass_scale*weighted_area;
227  ebuff_s[0][2] = mass_scale*weighted_area;
228  }
229  else if(nnodes==4){
230  ebuff[0][0] = interpolationCoefficients[0];
231  ebuff[0][1] = interpolationCoefficients[1];
232  ebuff[0][2] = interpolationCoefficients[2];
233  ebuff[0][3] = interpolationCoefficients[3];
234 
235  ebuff_counter[0][0] = interpolationCoefficients[0];
236  ebuff_counter[0][1] = interpolationCoefficients[1];
237  ebuff_counter[0][2] = interpolationCoefficients[2];
238  ebuff_counter[0][3] = interpolationCoefficients[3];
239 
240  ebuff_s[0][0] = mass_scale*weighted_area;
241  ebuff_s[0][1] = mass_scale*weighted_area;
242  ebuff_s[0][2] = mass_scale*weighted_area;
243  ebuff_s[0][3] = mass_scale*weighted_area;
244  }
245 }
246 
247 
279 template<class T, class S, class emi_index_rank>
280 inline void construct_direct_unique_both_operators(
281  SF::abstract_matrix<T, S>*& operator_unique_to_both_faces,
282  SF::abstract_matrix<T, S>*& operator_both_to_unique_face,
283  hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank, emi_index_rank>>& map_elem_uniqueFace_to_elem_bothface,
284  const hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank, emi_index_rank>>& map_elem_uniqueFace_to_elem_oneface,
285  const SF::vector<mesh_int_t>& vec_both_to_one_face,
286  const SF::meshdata<T, mesh_real_t>& emi_surfmesh_w_counter_face,
287  const SF::meshdata<T, mesh_real_t>& emi_surfmesh_unique_face,
288  int max_row_entries_emi,
289  int dpn,
290  typename SF::abstract_vector<T, S>::ltype alg_surface_type)
291 {
292  // Global/local dimensions for the direct unique-face <-> both-face operators.
293  T M = emi_surfmesh_w_counter_face.g_numelem;
294  T m = emi_surfmesh_w_counter_face.l_numelem;
295  T M_unique_face = emi_surfmesh_unique_face.g_numelem;
296  T m_unique_face = emi_surfmesh_unique_face.l_numelem;
297 
298  int rank = 0;
299  int comm_size = 0;
300  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
301  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
302 
303  SF::vector<long int> layout;
304  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout, emi_surfmesh_w_counter_face.comm);
305  T m_l = layout[rank];
306 
307  SF::vector<long int> layout_unique_face;
308  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_face, emi_surfmesh_unique_face.comm);
309  T m_unique_face_l = layout_unique_face[rank];
310 
311  // Rebuild the direct mapping from scratch every time the operators are assembled.
312  map_elem_uniqueFace_to_elem_bothface.clear();
313 
314  // Gather the local both->one lookup so every rank can recover the matching
315  // both-face indices for remote one-face owners.
316  std::vector<int> both_counts(comm_size, 0), both_displs(comm_size, 0);
317  int local_both_count = static_cast<int>(vec_both_to_one_face.size());
318  MPI_Allgather(&local_both_count, 1, MPI_INT, both_counts.data(), 1, MPI_INT,
319  emi_surfmesh_w_counter_face.comm);
320 
321  int total_both_count = 0;
322  for (int i = 0; i < comm_size; i++) {
323  both_displs[i] = total_both_count;
324  total_both_count += both_counts[i];
325  }
326 
327  std::vector<mesh_int_t> all_both_to_one(total_both_count);
328  // Exchange raw bytes so the gather stays correct even if mesh_int_t differs
329  // from int on a given build.
330  std::vector<int> both_byte_counts(comm_size, 0), both_byte_displs(comm_size, 0);
331  for (int i = 0; i < comm_size; i++) {
332  both_byte_counts[i] = both_counts[i] * static_cast<int>(sizeof(mesh_int_t));
333  both_byte_displs[i] = both_displs[i] * static_cast<int>(sizeof(mesh_int_t));
334  }
335  const int local_both_bytes = local_both_count * static_cast<int>(sizeof(mesh_int_t));
336  MPI_Allgatherv(reinterpret_cast<const unsigned char*>(vec_both_to_one_face.data()),
337  local_both_bytes, MPI_BYTE,
338  reinterpret_cast<unsigned char*>(all_both_to_one.data()),
339  both_byte_counts.data(), both_byte_displs.data(), MPI_BYTE,
340  emi_surfmesh_w_counter_face.comm);
341 
342  // Build per-rank one->both lookup tables. A one-face entry can map to up to
343  // two both-face entries, so we store first and second occurrences explicitly.
344  std::vector<std::vector<mesh_int_t>> one_to_both_first(comm_size);
345  std::vector<std::vector<mesh_int_t>> one_to_both_second(comm_size);
346  for (int r = 0; r < comm_size; r++) {
347  mesh_int_t max_one = -1;
348  for (int i = 0; i < both_counts[r]; i++) {
349  mesh_int_t one_idx = all_both_to_one[both_displs[r] + i];
350  if (one_idx > max_one) max_one = one_idx;
351  }
352  if (max_one < 0) continue;
353 
354  one_to_both_first[r].assign(max_one + 1, -1);
355  one_to_both_second[r].assign(max_one + 1, -1);
356  for (int i = 0; i < both_counts[r]; i++) {
357  mesh_int_t one_idx = all_both_to_one[both_displs[r] + i];
358  if (one_idx < 0) continue;
359  if (one_to_both_first[r][one_idx] < 0) one_to_both_first[r][one_idx] = i;
360  else one_to_both_second[r][one_idx] = i;
361  }
362  }
363 
364  // Convert the stored unique->one relation into a direct unique->both relation.
365  for (const auto& [unique_idx, one_face_pair] : map_elem_uniqueFace_to_elem_oneface) {
366  const auto& first_oneface = one_face_pair.first;
367  const auto& second_oneface = one_face_pair.second;
368 
369  auto map_one_to_both = [&](const emi_index_rank& one_face_idx, bool use_second) {
370  emi_index_rank both_face_idx;
371  both_face_idx.index = -1;
372  both_face_idx.rank = -1;
373 
374  if (one_face_idx.rank < 0 || one_face_idx.rank >= comm_size || one_face_idx.index < 0) {
375  return both_face_idx;
376  }
377  if (one_face_idx.index >= static_cast<int>(one_to_both_first[one_face_idx.rank].size())) {
378  return both_face_idx;
379  }
380 
381  mesh_int_t both_idx = one_to_both_first[one_face_idx.rank][one_face_idx.index];
382  if (use_second &&
383  one_face_idx.index < static_cast<int>(one_to_both_second[one_face_idx.rank].size()) &&
384  one_to_both_second[one_face_idx.rank][one_face_idx.index] >= 0) {
385  both_idx = one_to_both_second[one_face_idx.rank][one_face_idx.index];
386  }
387 
388  if (both_idx >= 0) {
389  both_face_idx.index = both_idx;
390  both_face_idx.rank = one_face_idx.rank;
391  }
392  return both_face_idx;
393  };
394 
395  const bool same_oneface =
396  (first_oneface.index >= 0 && second_oneface.index >= 0 &&
397  first_oneface.index == second_oneface.index &&
398  first_oneface.rank == second_oneface.rank);
399 
400  map_elem_uniqueFace_to_elem_bothface[unique_idx] = std::make_pair(
401  map_one_to_both(first_oneface, false),
402  map_one_to_both(second_oneface, same_oneface));
403  }
404 
405  // Assemble the expansion operator unique -> both. Each unique-face column
406  // writes into one or two both-face rows, depending on ownership/splitting.
407  SF::init_matrix(&operator_unique_to_both_faces);
408  operator_unique_to_both_faces->init(M, M_unique_face, m, m_unique_face, m_l, max_row_entries_emi*10);
409  operator_unique_to_both_faces->zero();
410 
411  {
412  // Broadcast owner-provided row locations so each rank can insert its owned rows.
413  const int M_unique_global = static_cast<int>(emi_surfmesh_unique_face.g_numelem);
414  std::vector<int> first_rank(M_unique_global, -1), first_idx(M_unique_global, -1);
415  std::vector<int> second_rank(M_unique_global, -1), second_idx(M_unique_global, -1);
416 
417  for (const auto& [local_unique_idx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
418  int global_unique = static_cast<int>(layout_unique_face[rank] + local_unique_idx);
419  if (global_unique < 0 || global_unique >= M_unique_global) continue;
420 
421  first_rank[global_unique] = static_cast<int>(both_pair.first.rank);
422  first_idx[global_unique] = static_cast<int>(both_pair.first.index);
423  second_rank[global_unique] = static_cast<int>(both_pair.second.rank);
424  second_idx[global_unique] = static_cast<int>(both_pair.second.index);
425  }
426 
427  MPI_Allreduce(MPI_IN_PLACE, first_rank.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
428  MPI_Allreduce(MPI_IN_PLACE, first_idx.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
429  MPI_Allreduce(MPI_IN_PLACE, second_rank.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
430  MPI_Allreduce(MPI_IN_PLACE, second_idx.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
431 
432  SF::vector<SF_int> row_idx(1), col_idx(1);
433  SF::dmat<SF_real> ebuff(1, 1);
434  ebuff.assign(1, 1, 1.0);
435 
436  for (int gid = 0; gid < M_unique_global; gid++) {
437  col_idx[0] = gid;
438  if (first_rank[gid] == rank && first_idx[gid] >= 0) {
439  row_idx[0] = static_cast<SF_int>(layout[rank] + first_idx[gid]);
440  operator_unique_to_both_faces->set_values(row_idx, col_idx, ebuff.data(), false);
441  }
442  if (second_rank[gid] == rank && second_idx[gid] >= 0) {
443  row_idx[0] = static_cast<SF_int>(layout[rank] + second_idx[gid]);
444  operator_unique_to_both_faces->set_values(row_idx, col_idx, ebuff.data(), false);
445  }
446  }
447  operator_unique_to_both_faces->finish_assembly();
448  }
449 
450  // Assemble the restriction/averaging operator both -> unique from the direct map.
451  SF::init_matrix(&operator_both_to_unique_face);
452  operator_both_to_unique_face->init(M_unique_face, M, m_unique_face, m, m_unique_face_l, max_row_entries_emi*10);
453  operator_both_to_unique_face->zero();
454  assemble_map_both_to_unique(*operator_both_to_unique_face,
455  map_elem_uniqueFace_to_elem_bothface,
456  emi_surfmesh_unique_face,
457  emi_surfmesh_w_counter_face);
458 
459 #ifdef EMI_DEBUG_MESH
460  {
461  // Validate mapping coverage and check whether the direct operators compose
462  // to the identity on a constant unique-face field.
463  mesh_int_t local_valid_first = 0, local_valid_second = 0;
464  mesh_int_t local_same_column = 0;
465  SF::vector<long int> layout_both_dbg;
466  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout_both_dbg, emi_surfmesh_w_counter_face.comm);
467  SF::vector<long int> layout_unique_dbg;
468  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_dbg, emi_surfmesh_unique_face.comm);
469 
470  for (const auto& [uidx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
471  if (both_pair.first.index >= 0) local_valid_first++;
472  if (both_pair.second.index >= 0) local_valid_second++;
473 
474  if (both_pair.first.index >= 0 && both_pair.second.index >= 0 &&
475  both_pair.first.rank >= 0 && both_pair.second.rank >= 0) {
476  mesh_int_t global_first = layout_both_dbg[both_pair.first.rank] + both_pair.first.index;
477  mesh_int_t global_second = layout_both_dbg[both_pair.second.rank] + both_pair.second.index;
478  if (global_first == global_second) local_same_column++;
479  }
480  }
481  mesh_int_t global_valid_first = 0, global_valid_second = 0;
482  mesh_int_t global_same_column = 0;
483  unsigned long long local_valid_first_ull = static_cast<unsigned long long>(local_valid_first);
484  unsigned long long local_valid_second_ull = static_cast<unsigned long long>(local_valid_second);
485  unsigned long long local_same_column_ull = static_cast<unsigned long long>(local_same_column);
486  unsigned long long global_valid_first_ull = 0;
487  unsigned long long global_valid_second_ull = 0;
488  unsigned long long global_same_column_ull = 0;
489  MPI_Allreduce(&local_valid_first_ull, &global_valid_first_ull, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, emi_surfmesh_w_counter_face.comm);
490  MPI_Allreduce(&local_valid_second_ull, &global_valid_second_ull, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, emi_surfmesh_w_counter_face.comm);
491  MPI_Allreduce(&local_same_column_ull, &global_same_column_ull, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, emi_surfmesh_w_counter_face.comm);
492  global_valid_first = static_cast<mesh_int_t>(global_valid_first_ull);
493  global_valid_second = static_cast<mesh_int_t>(global_valid_second_ull);
494  global_same_column = static_cast<mesh_int_t>(global_same_column_ull);
495 
496  int dbg_rank = -1;
497  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &dbg_rank);
498  if (dbg_rank == 0) {
499  log_msg(NULL, 0, 0, "DEBUG map_unique_to_both: valid_first=%d valid_second=%d (global M=%zu, M_unique=%zu)",
500  (int)global_valid_first, (int)global_valid_second,
501  (size_t)emi_surfmesh_w_counter_face.g_numelem,
502  (size_t)emi_surfmesh_unique_face.g_numelem);
503  log_msg(NULL, 0, 0, "DEBUG map_unique_to_both: same_column=%d", (int)global_same_column);
504  }
505 
506  {
507  std::vector<char> present(emi_surfmesh_unique_face.l_numelem, 0);
508  for (const auto& [uidx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
509  if (uidx >= 0 && uidx < (mesh_int_t)present.size()) present[uidx] = 1;
510  }
511  std::vector<int> missing;
512  for (mesh_int_t i = 0; i < (mesh_int_t)present.size(); i++) {
513  if (!present[i]) missing.push_back((int)i);
514  }
515  int local_missing = (int)missing.size();
516  int comm_size_dbg = 0;
517  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size_dbg);
518  std::vector<int> counts(comm_size_dbg, 0), displs(comm_size_dbg, 0);
519  MPI_Gather(&local_missing, 1, MPI_INT,
520  dbg_rank == 0 ? counts.data() : nullptr, 1, MPI_INT,
521  0, emi_surfmesh_w_counter_face.comm);
522  if (dbg_rank == 0) {
523  int total = 0;
524  for (int r = 0; r < comm_size_dbg; r++) {
525  displs[r] = total;
526  total += counts[r];
527  }
528  std::vector<int> all_missing(total, -1);
529  MPI_Gatherv(missing.data(), local_missing, MPI_INT,
530  all_missing.data(), counts.data(), displs.data(), MPI_INT,
531  0, emi_surfmesh_w_counter_face.comm);
532  int offset = 0;
533  bool any = false;
534  for (int r = 0; r < comm_size_dbg; r++) {
535  if (counts[r] == 0) {
536  offset += counts[r];
537  continue;
538  }
539  any = true;
540  log_msg(NULL, 0, 0, "DEBUG unique missing: rank=%d count=%d", r, counts[r]);
541  int to_print = counts[r] < 10 ? counts[r] : 10;
542  for (int i = 0; i < to_print; i++) {
543  log_msg(NULL, 0, 0, "DEBUG unique missing: rank=%d local_unique=%d", r, all_missing[offset + i]);
544  }
545  offset += counts[r];
546  }
547  if (!any) log_msg(NULL, 0, 0, "DEBUG unique missing: none");
548  } else {
549  MPI_Gatherv(missing.data(), local_missing, MPI_INT,
550  nullptr, nullptr, nullptr, MPI_INT,
551  0, emi_surfmesh_w_counter_face.comm);
552  }
553  }
554 
555  SF::abstract_vector<T, S>* dbg_unique = nullptr;
556  SF::abstract_vector<T, S>* dbg_both = nullptr;
557  SF::abstract_vector<T, S>* dbg_back = nullptr;
558  SF::init_vector(&dbg_unique, emi_surfmesh_unique_face, dpn, alg_surface_type);
559  SF::init_vector(&dbg_both, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
560  SF::init_vector(&dbg_back, emi_surfmesh_unique_face, dpn, alg_surface_type);
561  dbg_unique->set(1.0);
562  operator_unique_to_both_faces->mult(*dbg_unique, *dbg_both);
563  operator_both_to_unique_face->mult(*dbg_both, *dbg_back);
564 
565  SF_real local_max_err = 0.0;
566  mesh_int_t local_half_count = 0;
567  SF_real* p = dbg_back->ptr();
568  for (mesh_int_t i = 0; i < dbg_back->lsize(); i++) {
569  SF_real err = std::abs(p[i] - 1.0);
570  if (err > local_max_err) local_max_err = err;
571  if (std::abs(p[i] - 0.5) < 1e-12) local_half_count++;
572  }
573  dbg_back->release_ptr(p);
574  SF_real global_max_err = 0.0;
575  mesh_int_t global_half_count = 0;
576  MPI_Allreduce(&local_max_err, &global_max_err, 1, MPI_DOUBLE, MPI_MAX, emi_surfmesh_w_counter_face.comm);
577  unsigned long long local_half_count_ull = static_cast<unsigned long long>(local_half_count);
578  unsigned long long global_half_count_ull = 0;
579  MPI_Allreduce(&local_half_count_ull, &global_half_count_ull, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, emi_surfmesh_w_counter_face.comm);
580  global_half_count = static_cast<mesh_int_t>(global_half_count_ull);
581  if (dbg_rank == 0) {
582  log_msg(NULL, 0, 0, "DEBUG map_unique_to_both consistency: max_err=%.6e half_count=%d",
583  (double)global_max_err, (int)global_half_count);
584  }
585 
586  {
587  const int max_print = 5;
588  const int fields = 9;
589  std::vector<int> local_buf(max_print * fields, -1);
590  int filled = 0;
591 
592  for (const auto& [local_unique_idx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
593  if (filled >= max_print) break;
594  const auto& first_both = both_pair.first;
595  const auto& second_both = both_pair.second;
596 
597  const bool first_valid = (first_both.index >= 0 && first_both.rank >= 0);
598  const bool second_valid = (second_both.index >= 0 && second_both.rank >= 0);
599  int count = (first_valid ? 1 : 0) + (second_valid ? 1 : 0);
600  if (count != 1) continue;
601 
602  mesh_int_t global_unique = layout_unique_dbg[dbg_rank] + local_unique_idx;
603  mesh_int_t global_first = first_valid ? (layout_both_dbg[first_both.rank] + first_both.index) : -1;
604  mesh_int_t global_second = second_valid ? (layout_both_dbg[second_both.rank] + second_both.index) : -1;
605 
606  int base = filled * fields;
607  local_buf[base + 0] = dbg_rank;
608  local_buf[base + 1] = (int)local_unique_idx;
609  local_buf[base + 2] = (int)global_unique;
610  local_buf[base + 3] = (int)first_both.rank;
611  local_buf[base + 4] = (int)first_both.index;
612  local_buf[base + 5] = (int)global_first;
613  local_buf[base + 6] = (int)second_both.rank;
614  local_buf[base + 7] = (int)second_both.index;
615  local_buf[base + 8] = (int)global_second;
616  filled++;
617  }
618 
619  int comm_size = 0;
620  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
621  std::vector<int> all_buf;
622  if (dbg_rank == 0) all_buf.resize(comm_size * max_print * fields, -1);
623 
624  MPI_Gather(local_buf.data(), max_print * fields, MPI_INT,
625  dbg_rank == 0 ? all_buf.data() : nullptr, max_print * fields, MPI_INT,
626  0, emi_surfmesh_w_counter_face.comm);
627 
628  if (dbg_rank == 0) {
629  for (int r = 0; r < comm_size; r++) {
630  for (int i = 0; i < max_print; i++) {
631  int base = (r * max_print + i) * fields;
632  if (all_buf[base + 1] < 0) continue;
633  log_msg(NULL, 0, 0,
634  "DEBUG map_unique_to_both bad: rank=%d local_unique=%d global_unique=%d "
635  "first(rank=%d idx=%d glob=%d) second(rank=%d idx=%d glob=%d)",
636  all_buf[base + 0], all_buf[base + 1], all_buf[base + 2],
637  all_buf[base + 3], all_buf[base + 4], all_buf[base + 5],
638  all_buf[base + 6], all_buf[base + 7], all_buf[base + 8]);
639  }
640  }
641  }
642  }
643 
644  delete dbg_unique;
645  delete dbg_both;
646  delete dbg_back;
647  }
648 #endif
649 }
650 
697 template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
698 inline void assemble_restrict_operator( abstract_matrix<T,S> & B,
699  abstract_matrix<T,S> & Bi,
700  abstract_matrix<T,S> & BsM,
701  SF::vector<mesh_int_t> elemTag_surface_mesh,
702  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>,
703  std::pair<mesh_int_t,mesh_int_t>> map_vertex_tag_to_dof_petsc,
704  hashmap::unordered_map<tuple_key,
705  std::pair<tuple_value, tuple_value>> line_face,
706  hashmap::unordered_map<tri_key,
707  std::pair<tri_value, tri_value>> tri_face,
708  hashmap::unordered_map<quad_key,
709  std::pair<quad_value, quad_value>> quad_face,
710  const meshdata<mesh_int_t,mesh_real_t> & surface_mesh,
711  const meshdata<mesh_int_t,mesh_real_t> & emi_mesh,
712  S mass_scale)
713 {
714  // we want to make sure that the element integrator fits the matrix
715  T row_dpn = 1; T col_dpn = 1;
716 
717  const vector<SF_int> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);// remove it
718  const SF::vector<mesh_int_t> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
719  hashmap::unordered_map<T,T> g2l_emi;
720  hashmap::unordered_map<T,T> l2g_emi;
721  for(size_t i=0; i<rnod_emi.size(); i++){
722  g2l_emi[rnod_emi[i]] = i;
723  l2g_emi[i] = rnod_emi[i];
724  }
725 
726 
727  const SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
728  hashmap::unordered_map<T,T> g2l;
729  hashmap::unordered_map<T,T> l2g;
730 
731  for(size_t i=0; i<rnod.size(); i++){
732  g2l[rnod[i]] = i;
733  l2g[i] = rnod[i];
734  }
735 
736  // allocate row / col index buffers
737  vector<SF_int> row_idx(SF_MAX_ELEM_NODES * row_dpn), row_idx_counter(SF_MAX_ELEM_NODES * row_dpn), col_idx(SF_MAX_ELEM_NODES * col_dpn), col_idx_counter(SF_MAX_ELEM_NODES * col_dpn);
738 
739  const SF::vector<mesh_int_t> & emi_surfmesh_elem = surface_mesh.get_numbering(SF::NBR_ELEM_REF);
740 
741  // allocate elem index buffer
742  dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
743  dmat<SF_real> ebuff_s(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
744  dmat<SF_real> ebuff_counter(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
745 
746  const T* con_emi = emi_mesh.con.data();
747  const T* con = surface_mesh.con.data();
748 
749  // start with assembly
750  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
751  {
752  T tag = surface_mesh.tag[eidx];
753 
754  // get the counter part
755  std::vector<int> elem_nodes;
756  std::vector<int> elem_nodes_old;
757  std::vector<int> petsc_first;
758  std::vector<int> petsc_second;
759 
760  T tag_first = 0;
761  T tag_second = 0;
762 
763  T rank_first = 0;
764  T rank_second = 0;
765 
766  T mem_first = 0;
767  T mem_second = 0;
768 
769  T sign = 1.0;
770 
771  for (int n = surface_mesh.dsp[eidx]; n < surface_mesh.dsp[eidx+1];n++)
772  {
773  T l_idx = surface_mesh.con[n];
774 
775  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
776  Index_tag_old = std::make_pair(l2g[l_idx],tag);
777  mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
778  elem_nodes.push_back(Index_new);
779  elem_nodes_old.push_back(l2g[l_idx]);
780  }
781 
782  std::sort(elem_nodes.begin(),elem_nodes.end());
783  if(elem_nodes.size()==2){
784  tuple_key key;
785 
786  key.v1 = elem_nodes[0];
787  key.v2 = elem_nodes[1];
788  std::pair<tuple_value, tuple_value> value = line_face[key];
789 
790  tag_first = value.first.tag;
791  tag_second = value.second.tag;
792 
793  rank_first = value.first.rank;
794  rank_second = value.second.rank;
795 
796  mem_first = value.first.mem;
797  mem_second = value.second.mem;
798  }
799  else if(elem_nodes.size()==3){
800  tri_key key;
801  key.v1 = elem_nodes[0];
802  key.v2 = elem_nodes[1];
803  key.v3 = elem_nodes[2];
804  std::pair<tri_value, tri_value> value = tri_face[key];
805 
806  tag_first = value.first.tag;
807  tag_second = value.second.tag;
808 
809  rank_first = value.first.rank;
810  rank_second = value.second.rank;
811 
812  mem_first = value.first.mem;
813  mem_second = value.second.mem;
814  }
815  else if(elem_nodes.size()==4){
816  quad_key key;
817  key.v1 = elem_nodes[0];
818  key.v2 = elem_nodes[1];
819  key.v3 = elem_nodes[2];
820  key.v4 = elem_nodes[3];
821  std::pair<quad_value, quad_value> value = quad_face[key];
822 
823  tag_first = value.first.tag;
824  tag_second = value.second.tag;
825 
826  rank_first = value.first.rank;
827  rank_second = value.second.rank;
828 
829  mem_first = value.first.mem;
830  mem_second = value.second.mem;
831  }
832 
833  if(tag != tag_first) std::swap(tag_first, tag_second);
834 
835  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
836 
837  for (int indx = 0; indx < elem_nodes_old.size(); ++indx)
838  {
839  std::pair <mesh_int_t,mesh_int_t> Index_tag_old_first;
840  Index_tag_old_first = std::make_pair(elem_nodes_old[indx],tag_first);
841 
842  auto it_first = map_vertex_tag_to_dof_petsc.find(Index_tag_old_first);
843  if (it_first == map_vertex_tag_to_dof_petsc.end()) {
844  std::cerr << "ERROR: tag_first=" << tag_first << " not found for vertex=" << elem_nodes_old[indx] << std::endl;
845  throw std::runtime_error("PETSc index not found for first tag");
846  }
847  std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_first = it_first->second;
848  mesh_int_t newIndex_first = newIndex_petsc_first.first;
849  petsc_first.push_back(newIndex_petsc_first.second);
850 
851  std::pair <mesh_int_t,mesh_int_t> Index_tag_old_second;
852  Index_tag_old_second = std::make_pair(elem_nodes_old[indx],tag_second);
853 
854  auto it_second = map_vertex_tag_to_dof_petsc.find(Index_tag_old_second);
855  if (it_second == map_vertex_tag_to_dof_petsc.end()) {
856  std::cerr << "ERROR: tag_second=" << tag_second << " not found for vertex=" << elem_nodes_old[indx] << std::endl;
857  throw std::runtime_error("PETSc index not found for counter-face tag");
858  }
859  std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_second = it_second->second;
860  mesh_int_t newIndex_second = newIndex_petsc_second.first;
861  petsc_second.push_back(newIndex_petsc_second.second);
862 
863  // elemTag_surface_mesh[eidx] == 1 means that the tag belongs to the extracellular region
864  if(mem_first==1 and elemTag_surface_mesh[eidx]==1) // condition for the membrane which the sign is always negative on extracellular side
865  sign = -1.0;
866  else if (mem_first==2 and tag<tag_second) // condition for the gap juntion which the sign is always negative for smaller tag
867  sign = -1.0;
868  }
869 
870  SF_int nnodes = elem_nodes.size();
871  vector<SF::Point> face_coordinates(nnodes);
872  // calculate row/col indices of entries
873  row_idx.resize(nnodes*row_dpn);
874  row_idx_counter.resize(nnodes*row_dpn);
875  col_idx.resize(nnodes*col_dpn);
876  col_idx_counter.resize(nnodes*col_dpn);
877 
878  {
879  for(SF_int i=0; i<nnodes; i++){
880  // all the row sets with element index
881  for(short j=0; j<row_dpn; j++){
882  row_idx.data()[i*row_dpn + j] = emi_surfmesh_elem[eidx]*row_dpn;
883  }
884 
885  for(short j=0; j<col_dpn; j++){
886  col_idx[i*col_dpn + j] = (petsc_first[i])*col_dpn + j;
887  col_idx_counter[i*col_dpn + j] = (petsc_second[i])*col_dpn + j;
888  }
889  }
890  }
891 
892  // assign ebuff and ebuff_s
893  ebuff.assign(nnodes, nnodes, 0.0);
894  ebuff_s.assign(nnodes, nnodes, 0.0);
895 
896  ebuff_counter.assign(nnodes, nnodes, 0.0);
897 
898  for (int n = surface_mesh.dsp[eidx], i = 0; n < surface_mesh.dsp[eidx+1];n++,i++)
899  {
900  T l_idx = surface_mesh.con[n];
901 
902  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
903  Index_tag_old = std::make_pair(l2g[l_idx],tag);
904  mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
905 
906  double x = emi_mesh.xyz[g2l_emi[Index_new]*3+0];
907  double y = emi_mesh.xyz[g2l_emi[Index_new]*3+1];
908  double z = emi_mesh.xyz[g2l_emi[Index_new]*3+2];
909 
910  face_coordinates[i].x = x;
911  face_coordinates[i].y = y;
912  face_coordinates[i].z = z;
913  }
914  compute_integrate_matrix_barycentric(face_coordinates, nnodes, ebuff, ebuff_s, ebuff_counter, mass_scale);
915 
916  // add values into system matrix
917  const bool add = true;
918 
919  // we simply build Bi in the same way as B, but without the sign change, so we end up with Bi*Iij_stim = Ib_stim = Ii + Ij
920  Bi.set_values(row_idx, col_idx, ebuff.data(), add);
921  Bi.set_values(row_idx, col_idx_counter, ebuff_counter.data(), add);
922  // the operator matrix in order to project potential of dofs to the barycenteric of each face as transmembrane (vm = ui - uj)
923  ebuff*=(sign);
924  B.set_values(row_idx, col_idx, ebuff.data(), add);
925  ebuff_counter*=(-sign);
926  B.set_values(row_idx, col_idx_counter, ebuff_counter.data(), add); // counter face
927 
928  // the operator matrix in order to project the current density to the dofs of the face
929  ebuff_s*=(sign);
930  BsM.set_values(col_idx, row_idx, ebuff_s.data(), add);
931  }
932  // finish assembly and progress output
933  B.finish_assembly();
934  Bi.finish_assembly();
935  BsM.finish_assembly();
936 }
937 
954 template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
955 inline void assemble_lhs_emi(abstract_matrix<T,S> & mat,
956  abstract_matrix<T,S> & mat_surf,
957  const meshdata<mesh_int_t,mesh_real_t> & emi_mesh,
958  const meshdata<mesh_int_t,mesh_real_t> & surface_mesh,
959  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>> map_vertex_tag_to_dof_petsc,
960  hashmap::unordered_map<tuple_key,
961  std::pair<tuple_value, tuple_value>> line_face,
962  hashmap::unordered_map<tri_key,
963  std::pair<tri_value, tri_value>> tri_face,
964  hashmap::unordered_map<quad_key,
965  std::pair<quad_value, quad_value>> quad_face,
966  matrix_integrator<mesh_int_t,mesh_real_t> & stiffness_integrator,
967  matrix_integrator<mesh_int_t, mesh_real_t> & mass_integrator,
968  S stiffness_scale,
969  S mass_scale)
970 {
971  // we want to make sure that the element integrator fits the matrix
972  T row_dpn = 1; T col_dpn = 1;
973 
974  // allocate row / col index buffers
975  vector<SF_int> idx(SF_MAX_ELEM_NODES * col_dpn), idx_counter(SF_MAX_ELEM_NODES * col_dpn);
976  vector<T> row_idx(SF_MAX_ELEM_NODES * row_dpn), col_idx(SF_MAX_ELEM_NODES * col_dpn);
977 
978  // allocate elem index buffer
979  dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * col_dpn, SF_MAX_ELEM_NODES * col_dpn);
980  dmat<SF_real> ebuff_mass(SF_MAX_ELEM_NODES * col_dpn, SF_MAX_ELEM_NODES * col_dpn);
981  // Assemble stiffness matrix (K)
982  const vector<mesh_int_t> & stiffness_petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);
983  element_view<mesh_int_t, mesh_real_t> stiffness_view(emi_mesh, NBR_PETSC);
984 
985  // Assembly of global stiffness matrix into lhs matrix
986  for(size_t eidx=0; eidx < emi_mesh.l_numelem; eidx++)
987  {
988  // set element view to current element
989  stiffness_view.set_elem(eidx);
990 
991  // calculate row/col indices of entries
992  mesh_int_t nnodes = stiffness_view.num_nodes();
993 
994  row_idx.resize(nnodes*row_dpn);
995  col_idx.resize(nnodes*col_dpn);
996  canonic_indices<mesh_int_t,SF_int>(stiffness_view.nodes(), stiffness_petsc_nbr.data(), nnodes, row_dpn, row_idx.data());
997  canonic_indices<mesh_int_t,SF_int>(stiffness_view.nodes(), stiffness_petsc_nbr.data(), nnodes, col_dpn, col_idx.data());
998 
999  // call integrator
1000  stiffness_integrator(stiffness_view, ebuff);
1001  ebuff *= stiffness_scale;
1002 
1003  // add values into system matrix
1004  const bool add = true;
1005  mat.set_values(row_idx, col_idx, ebuff.data(), add);
1006  }
1007 
1008  // Assemble the global surface mass matrix (M).
1009  // Each local mass matrix computed for a surface element is mapped into the global system
1010  // using the corresponding global row and column indices of that element’s nodes.
1011 
1012  const vector<SF_int> & ref_nbr = surface_mesh.get_numbering(SF::NBR_REF);
1013  const vector<SF_int> & petsc_nbr = surface_mesh.get_numbering(SF::NBR_PETSC);
1014 
1015  const SF::vector<mesh_int_t> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
1016  hashmap::unordered_map<T,T> g2l_emi;
1017  hashmap::unordered_map<T,T> l2g_emi;
1018  for(size_t i=0; i<rnod_emi.size(); i++){
1019  g2l_emi[rnod_emi[i]] = i;
1020  l2g_emi[i] = rnod_emi[i];
1021  }
1022  const SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
1023  hashmap::unordered_map<T,T> g2l;
1024  hashmap::unordered_map<T,T> l2g;
1025 
1026  for(size_t i=0; i<rnod.size(); i++){
1027  g2l[rnod[i]] = i;
1028  l2g[i] = rnod[i];
1029  }
1030 
1031  // start with assembly mass matrix into lhs matrix
1032  element_view<mesh_int_t, mesh_real_t> view(surface_mesh, NBR_PETSC);
1033  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
1034  {
1035  // set element view to current element
1036  view.set_elem(eidx);
1037  mesh_int_t tag = surface_mesh.tag[eidx]; // selected face
1038 
1039  // get the counter part
1040  std::vector<int> elem_nodes;
1041  std::vector<int> elem_nodes_old;
1042  std::vector<int> elem_nodes_first;
1043  std::vector<int> elem_nodes_second;
1044 
1045  std::vector<int> petsc_first;
1046  std::vector<int> petsc_second;
1047 
1048  T tag_first = 0;
1049  T tag_second = 0;
1050 
1051  for (int n = surface_mesh.dsp[eidx]; n < surface_mesh.dsp[eidx+1];n++)
1052  {
1053  T l_idx = surface_mesh.con[n];
1054 
1055  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
1056  Index_tag_old = std::make_pair(l2g[l_idx],tag);
1057  mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
1058  elem_nodes.push_back(Index_new);
1059  elem_nodes_old.push_back(l2g[l_idx]);
1060  }
1061 
1062  std::sort(elem_nodes.begin(),elem_nodes.end());
1063  if(elem_nodes.size()==2){
1064  tuple_key key;
1065 
1066  key.v1 = elem_nodes[0];
1067  key.v2 = elem_nodes[1];
1068  std::pair<tuple_value, tuple_value> value = line_face[key];
1069 
1070  tag_first = value.first.tag;
1071  tag_second = value.second.tag;
1072  }
1073  else if(elem_nodes.size()==3){
1074  tri_key key;
1075  key.v1 = elem_nodes[0];
1076  key.v2 = elem_nodes[1];
1077  key.v3 = elem_nodes[2];
1078  std::pair<tri_value, tri_value> value = tri_face[key];
1079 
1080  tag_first = value.first.tag;
1081  tag_second = value.second.tag;
1082  }
1083  else if(elem_nodes.size()==4){
1084  quad_key key;
1085  key.v1 = elem_nodes[0];
1086  key.v2 = elem_nodes[1];
1087  key.v3 = elem_nodes[2];
1088  key.v4 = elem_nodes[3];
1089  std::pair<quad_value, quad_value> value = quad_face[key];
1090 
1091  tag_first = value.first.tag;
1092  tag_second = value.second.tag;
1093  }
1094 
1095  if(tag != tag_first) std::swap(tag_first, tag_second);
1096 
1097  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
1098 
1099  for (int indx = 0; indx < elem_nodes_old.size(); ++indx)
1100  {
1101  std::pair <mesh_int_t,mesh_int_t> Index_tag_old_first;
1102  Index_tag_old_first = std::make_pair(elem_nodes_old[indx],tag_first);
1103  std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_first = map_vertex_tag_to_dof_petsc[Index_tag_old_first];
1104  mesh_int_t newIndex_first = newIndex_petsc_first.first;
1105  petsc_first.push_back(newIndex_petsc_first.second);
1106 
1107  std::pair <mesh_int_t,mesh_int_t> Index_tag_old_second;
1108  Index_tag_old_second = std::make_pair(elem_nodes_old[indx],tag_second);
1109  std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_second = map_vertex_tag_to_dof_petsc[Index_tag_old_second];
1110  mesh_int_t newIndex_second = newIndex_petsc_second.first;
1111  petsc_second.push_back(newIndex_petsc_second.second);
1112  }
1113 
1114  // calculate row/col indices of entries
1115  SF_int nnodes = view.num_nodes();
1116  idx.resize(nnodes*row_dpn);
1117  idx_counter.resize(nnodes*col_dpn);
1118 
1119  {
1120  for(T i=0; i<nnodes; i++)
1121  {
1122  for(short j=0; j<row_dpn; j++){
1123  idx[i*row_dpn + j] = (petsc_first[i])*row_dpn + j;
1124  }
1125  for(short j=0; j<col_dpn; j++){
1126  idx_counter[i*col_dpn + j] = (petsc_second[i])*col_dpn + j;
1127  }
1128  }
1129  }
1130 
1131  // add values into system matrix
1132  const bool add = true;
1133 
1134  // call integrator
1135  mass_integrator(view, ebuff);
1136  ebuff_mass.assign(nnodes, nnodes, 0.0);
1137  ebuff_mass = ebuff;
1138  // Note: surface mesh by default has both face and counter face;
1139  // to avoid the dublication, factor 0.5 is necessary
1140  ebuff*=0.5 * mass_scale;
1141  ebuff_mass*=0.5;
1142 
1143  mat.set_values(idx, idx, ebuff.data(), add);
1144  mat.set_values(idx_counter, idx_counter, ebuff.data(), add);
1145 
1146  mat_surf.set_values(idx, idx, ebuff_mass.data(), add);
1147  mat_surf.set_values(idx_counter, idx_counter, ebuff_mass.data(), add);
1148 
1149  ebuff*=(-1.0);
1150  ebuff_mass*=(-1.0);
1151 
1152  mat.set_values(idx, idx_counter, ebuff.data(), add);
1153  mat.set_values(idx_counter, idx, ebuff.data(), add);
1154 
1155  mat_surf.set_values(idx, idx_counter, ebuff_mass.data(), add);
1156  mat_surf.set_values(idx_counter, idx, ebuff_mass.data(), add);
1157 
1158  // Each face contributes the following 2x2 block structure:
1159  //
1160  // s* [ + M - M ;
1161  // - M + M ]
1162  //
1163  // with s = 0.5 * alpha,
1164  // where M is the local face mass matrix.
1165 
1166  }
1167  // finish assembly and progress output
1168  mat.finish_assembly();
1169  mat_surf.finish_assembly();
1170 }
1171 
1189 template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
1190 inline void assign_resting_potential_from_ionic_models_on_myocyte(abstract_matrix<T,S> & Vm_emi,
1191  abstract_vector<T, S>* vb,
1192  SF::vector<mesh_int_t> elemTag_emi_mesh,
1193  hashmap::unordered_map<std::pair<T,T>,
1194  std::pair<T,T>> map_vertex_tag_to_dof_petsc,
1195  hashmap::unordered_map<tuple_key,
1196  std::pair<tuple_value, tuple_value>> line_face,
1197  hashmap::unordered_map<tri_key,
1198  std::pair<tri_value, tri_value>> tri_face,
1199  hashmap::unordered_map<quad_key,
1200  std::pair<quad_value, quad_value>> quad_face,
1201  const meshdata<T,mesh_real_t> & surface_mesh,
1202  const meshdata<T,mesh_real_t> & emi_mesh)
1203 {
1204  // we want to make sure that the element integrator fits the matrix
1205  T row_dpn = 1; T col_dpn = 1;
1206 
1207  const vector<SF_int> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);// remove it
1208  const SF::vector<T> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
1209  hashmap::unordered_map<T,T> g2l_emi;
1210  hashmap::unordered_map<T,T> l2g_emi;
1211  for(size_t i=0; i<rnod_emi.size(); i++){
1212  g2l_emi[rnod_emi[i]] = i;
1213  l2g_emi[i] = rnod_emi[i];
1214  }
1215 
1216  const SF::vector<T> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
1217  const vector<SF_int> & petsc_nbr_surface = surface_mesh.get_numbering(NBR_PETSC);// remove it
1218 
1219  hashmap::unordered_map<T,T> g2l;
1220  hashmap::unordered_map<T,T> l2g;
1221 
1222  hashmap::unordered_map<T,SF_real> tag_2_vm;
1223 
1224  for(size_t i=0; i<rnod.size(); i++){
1225  g2l[rnod[i]] = i;
1226  l2g[i] = rnod[i];
1227  }
1228  auto vb_data = vb->const_ptr();
1229 
1230  // start with assembly
1231  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
1232  {
1233  T tag = surface_mesh.tag[eidx];
1234 
1235  // get the counter part
1236  std::vector<int> elem_nodes;
1237  bool tag_current;
1238  T tag_first = 0;
1239  T tag_second = 0;
1240  T mem_first = 0;
1241  T mem_second = 0;
1242 
1243  T sign = 1.0;
1244 
1245  for (int n = surface_mesh.dsp[eidx], i = 0; n < surface_mesh.dsp[eidx+1];n++,i++)
1246  {
1247  T l_idx = surface_mesh.con[n];
1248  std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
1249  Index_tag_old = std::make_pair(l2g[l_idx],tag);
1250  mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
1251  elem_nodes.push_back(Index_new);
1252  }
1253 
1254  std::sort(elem_nodes.begin(),elem_nodes.end());
1255  if(elem_nodes.size()==2){
1256  tuple_key key;
1257 
1258  key.v1 = elem_nodes[0];
1259  key.v2 = elem_nodes[1];
1260  std::pair<tuple_value, tuple_value> value = line_face[key];
1261 
1262  tag_first = value.first.tag;
1263  tag_second = value.second.tag;
1264 
1265  mem_first = value.first.mem;
1266  mem_second = value.second.mem;
1267  }
1268  else if(elem_nodes.size()==3){
1269  tri_key key;
1270  key.v1 = elem_nodes[0];
1271  key.v2 = elem_nodes[1];
1272  key.v3 = elem_nodes[2];
1273  std::pair<tri_value, tri_value> value = tri_face[key];
1274 
1275  tag_first = value.first.tag;
1276  tag_second = value.second.tag;
1277 
1278  mem_first = value.first.mem;
1279  mem_second = value.second.mem;
1280  }
1281  else if(elem_nodes.size()==4){
1282  quad_key key;
1283  key.v1 = elem_nodes[0];
1284  key.v2 = elem_nodes[1];
1285  key.v3 = elem_nodes[2];
1286  key.v4 = elem_nodes[3];
1287  std::pair<quad_value, quad_value> value = quad_face[key];
1288 
1289  tag_first = value.first.tag;
1290  tag_second = value.second.tag;
1291 
1292  mem_first = value.first.mem;
1293  mem_second = value.second.mem;
1294  }
1295 
1296  SF_int nnodes = elem_nodes.size();
1297  for (int i = 0; i < nnodes; ++i)
1298  {
1299  if(mem_first==1 || mem_second==1){
1300  tag_2_vm[tag_first] = vb_data[eidx];
1301  tag_2_vm[tag_second] = vb_data[eidx];
1302  }
1303  }
1304  }
1305 
1306  vb->const_release_ptr(vb_data);
1307 
1308  // assign all the tags
1309  {
1310  // we want to make sure that the element integrator fits the matrix
1311  mesh_int_t row_dpn = 1, col_dpn = 1;
1312 
1313  // allocate row and value
1314  vector<T> row_idx(SF_MAX_ELEM_NODES * row_dpn);
1315  vector<T> col_idx(SF_MAX_ELEM_NODES * row_dpn);
1316  vector<SF_real> val_idx(SF_MAX_ELEM_NODES * row_dpn);
1317 
1318  // allocate elem index buffer
1319  dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
1320 
1321  const vector<mesh_int_t> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);
1322 
1323  // start with assembly
1324  element_view<mesh_int_t, mesh_real_t> view(emi_mesh, NBR_PETSC);
1325  for(size_t eidx=0; eidx < emi_mesh.l_numelem; eidx++)
1326  {
1327  // set element view to current element
1328  view.set_elem(eidx);
1329  T tag = emi_mesh.tag[eidx];
1330  // calculate row/col indices of entries
1331  mesh_int_t nnodes = view.num_nodes();
1332 
1333  row_idx.resize(nnodes*row_dpn);val_idx.resize(nnodes*row_dpn);
1334  col_idx.resize(nnodes*row_dpn);
1335  ebuff.assign(nnodes, nnodes, 0.0);
1336  canonic_indices<mesh_int_t,SF_int>(view.nodes(), petsc_nbr.data(), nnodes, row_dpn, row_idx.data());
1337  canonic_indices<mesh_int_t,SF_int>(view.nodes(), petsc_nbr.data(), nnodes, col_dpn, col_idx.data());
1338 
1339  for (int i = 0; i < nnodes; ++i)
1340  {
1341  ebuff[i][i] = 0;
1342  // elemTag_emi_mesh[eidx]==2 means that the tag belongs to the intracellular region
1343  if(elemTag_emi_mesh[eidx]==2){
1344  ebuff[i][i] = tag_2_vm[tag];
1345  }
1346  }
1347  Vm_emi.set_values(row_idx, col_idx, ebuff.data(), false);
1348  }
1349  Vm_emi.finish_assembly();
1350  }
1351 }
1352 
1353 }
1354 
1355 #endif
1356 #endif
sign
#define sign(x)
Definition: ION_IF.h:57
SF_container.h
Basic containers.
mesh_int_t
int mesh_int_t
Definition: SF_container.h:46
SF_fem_utils.h
FEM utilities.
SF_MAX_ELEM_NODES
#define SF_MAX_ELEM_NODES
max #nodes defining an element
Definition: SF_fem_utils.h:40
SF_real
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
SF_int
std::int32_t SF_int
Use the general std::int32_t as int type.
Definition: SF_globals.h:37
SF_init.h
SF::abstract_matrix
Definition: SF_abstract_matrix.h:45
SF::abstract_matrix::mult
virtual void mult(const abstract_vector< T, S > &x, abstract_vector< T, S > &b) const =0
SF::abstract_matrix::finish_assembly
virtual void finish_assembly()=0
SF::abstract_matrix::zero
virtual void zero()=0
SF::abstract_matrix::init
virtual void init(T iNRows, T iNCols, T ilrows, T ilcols, T loc_offset, T mxent)
Definition: SF_abstract_matrix.h:87
SF::abstract_matrix::set_values
virtual void set_values(const vector< T > &row_idx, const vector< T > &col_idx, const vector< S > &vals, bool add)=0
SF::abstract_vector
Definition: SF_abstract_vector.h:54
SF::abstract_vector::ptr
virtual S * ptr()=0
SF::abstract_vector::release_ptr
virtual void release_ptr(S *&p)=0
SF::abstract_vector::lsize
virtual T lsize() const =0
SF::abstract_vector::set
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false, const bool local=false)=0
SF::abstract_vector::ltype
ltype
Definition: SF_abstract_vector.h:59
SF::dmat
Dense matrix class.
Definition: dense_mat.hpp:43
SF::meshdata
The mesh storage class. It contains both element and vertex data.
Definition: SF_container.h:396
SF::meshdata::l_numelem
size_t l_numelem
local number of elements
Definition: SF_container.h:399
SF::meshdata::g_numelem
size_t g_numelem
global number of elements
Definition: SF_container.h:398
SF::meshdata::comm
MPI_Comm comm
the parallel mesh is defined on a MPI world
Definition: SF_container.h:404
SF::vector::size
size_t size() const
The current size of the vector.
Definition: SF_vector.h:104
SF::vector::assign
void assign(InputIterator s, InputIterator e)
Assign a memory range.
Definition: SF_vector.h:161
SF::vector::data
T * data()
Pointer to the vector's start.
Definition: SF_vector.h:91
hashmap::unordered_map
Definition: hashmap.hpp:253
log_msg
#define log_msg(F, L, O,...)
Definition: filament.h:8
SF
Definition: dense_mat.hpp:34
SF::mag
double mag(const Point &vect)
vector magnitude
Definition: SF_container.h:111
SF::inner_prod
double inner_prod(const Point &a, const Point &b)
Definition: SF_container.h:90
SF::count
void count(const vector< T > &data, vector< S > &cnt)
Count number of occurrences of indices.
Definition: SF_vector.h:332
SF::init_vector
void init_vector(SF::abstract_vector< T, S > **vec)
Definition: SF_init.h:99
SF::init_matrix
void init_matrix(SF::abstract_matrix< T, S > **mat)
Definition: SF_init.h:199
SF::cross
Point cross(const Point &a, const Point &b)
cross product
Definition: SF_container.h:84
SF::clamp
V clamp(const V val, const W start, const W end)
Clamp a value into an interval [start, end].
Definition: SF_base.h:63
SF::NBR_PETSC
@ NBR_PETSC
PETSc numbering of nodes.
Definition: SF_container.h:203
SF::NBR_ELEM_REF
@ NBR_ELEM_REF
The element numbering of the reference mesh (the one stored on HD).
Definition: SF_container.h:204
SF::NBR_REF
@ NBR_REF
The nodal numbering of the reference mesh (the one stored on HD).
Definition: SF_container.h:201
opencarp::Point
vec3< POINT_REAL > Point
Definition: vect.h:93
SF::Point
Point and vector struct.
Definition: SF_container.h:65
SF::Point::y
double y
Definition: SF_container.h:67
SF::Point::z
double z
Definition: SF_container.h:68
SF::Point::x
double x
Definition: SF_container.h:66
Generated by doxygen 1.9.1