openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
SF_mesh_io_emi.h
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
2 // SPDX-License-Identifier: Apache-2.0
3 
11 #if WITH_EMI_MODEL
12 
13 #ifndef _SF_MESH_IO_EMI
14 #define _SF_MESH_IO_EMI
15 
16 #include <cmath>
17 
18 #include "SF_mesh_io.h"
19 #include "backend.h"
20 
21 namespace SF {
22 
38 template<class T, class S>
39 inline void insert_points_ptsData_to_dof( const meshdata<T, S> & mesh_original, meshdata<T, S> & mesh,
40  const SF_nbr numbering,
41  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
42  hashmap::unordered_map<T,T> & vertex2ptsdata,
43  hashmap::unordered_map<T,T> & dof2ptsData)
44 {
45  MPI_Comm comm = mesh_original.comm;
46 
47  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
48  hashmap::unordered_map<T,T> g2l_orig;
49  g2l_orig.reserve(rnod_original.size());
50  for(size_t i=0; i<rnod_original.size(); i++){
51  g2l_orig[rnod_original[i]] = i;
52  }
53 
54  const SF::vector<mesh_int_t> & rnod = mesh.get_numbering(numbering);
55  hashmap::unordered_map<T,T> g2l;
56  hashmap::unordered_map<T,T> l2g;
57  g2l.reserve(rnod.size());
58  l2g.reserve(rnod.size());
59  for(size_t i=0; i<rnod.size(); i++){
60  g2l[rnod[i]] = i;
61  l2g[i] = rnod[i];
62  }
63 
64  mesh_int_t gmax = global_max(rnod, comm);
65  mesh.g_numpts = gmax+1;
66 
67  mesh.xyz.resize(mesh.l_numpts*3);
68 
69  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
70  T tag = mesh.tag[eidx];
71 
72  for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
73  {
74  T g = l2g[mesh.con[n]];
75  T old_g = dof2vertex[g];
76 
77  const auto local_idx = g2l[g];
78  const auto orig_idx = g2l_orig[old_g];
79 
80  mesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
81  mesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
82  mesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
83  dof2ptsData[g] = vertex2ptsdata[old_g];
84 
85  if (!std::isfinite(mesh.xyz[local_idx*3+0]) ||
86  !std::isfinite(mesh.xyz[local_idx*3+1]) ||
87  !std::isfinite(mesh.xyz[local_idx*3+2]))
88  {
89  printf("NaN detected in coordinates of the EMI mesh after decoupling interfaces on EMI mesh while recording pts data on the EMI mesh, on Tag=%lld\n",
90  static_cast<long long>(tag));
91  EXIT(EXIT_FAILURE);
92  }
93  }
94  }
95 }
96 
112 template<class T, class S>
113 inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh,
114  const SF_nbr numbering,
115  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
116  hashmap::unordered_set<int> & extra_tags,
117  SF::vector<mesh_int_t> & elemTag_surface)
118 {
119  MPI_Comm comm = mesh_original.comm;
120 
121  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
122  hashmap::unordered_map<T,T> g2l_orig;
123  g2l_orig.reserve(rnod_original.size());
124  for(size_t i=0; i<rnod_original.size(); i++){
125  g2l_orig[rnod_original[i]] = i;
126  }
127 
128  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
129 
130  hashmap::unordered_map<T,T> g2l_surf;
131  hashmap::unordered_map<T,T> l2g_surf;
132  g2l_surf.reserve(rnod_surf.size());
133  l2g_surf.reserve(rnod_surf.size());
134  for(size_t i=0; i<rnod_surf.size(); i++){
135  g2l_surf[rnod_surf[i]] = i;
136  l2g_surf[i] = rnod_surf[i];
137  }
138 
139  mesh_int_t gmax = global_max(rnod_surf, comm);
140  surfmesh.g_numpts = gmax+1;
141 
142  surfmesh.xyz.resize(surfmesh.l_numpts*3);
143 
144  elemTag_surface.resize(surfmesh.l_numelem);
145  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
146  T tag = surfmesh.tag[eidx];
147  elemTag_surface[eidx] = 1; // assigned 1 for for extracellular located on extracellular side
148  if(extra_tags.find(tag) == extra_tags.end())
149  elemTag_surface[eidx] = 2; // assigned 2 for for intracellular located on myocyte
150 
151  for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
152  {
153  T g = l2g_surf[surfmesh.con[n]];
154  T old_g = dof2vertex[g];
155 
156  const auto local_idx = g2l_surf[g];
157  const auto orig_idx = g2l_orig[old_g];
158 
159  surfmesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
160  surfmesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
161  surfmesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
162 
163  if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) ||
164  !std::isfinite(surfmesh.xyz[local_idx*3+1]) ||
165  !std::isfinite(surfmesh.xyz[local_idx*3+2]))
166  {
167  fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%lld \n",
168  static_cast<long long>(tag));
169  exit(1);
170  }
171  }
172  }
173 }
174 
186 template<class T, class S>
187 inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh,
188  const SF_nbr numbering,
189  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex)
190 {
191  MPI_Comm comm = mesh_original.comm;
192 
193  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
194  hashmap::unordered_map<T,T> g2l_orig;
195  g2l_orig.reserve(rnod_original.size());
196  for(size_t i=0; i<rnod_original.size(); i++){
197  g2l_orig[rnod_original[i]] = i;
198  }
199 
200  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
201 
202  hashmap::unordered_map<T,T> g2l_surf;
203  hashmap::unordered_map<T,T> l2g_surf;
204  g2l_surf.reserve(rnod_surf.size());
205  l2g_surf.reserve(rnod_surf.size());
206  for(size_t i=0; i<rnod_surf.size(); i++){
207  g2l_surf[rnod_surf[i]] = i;
208  l2g_surf[i] = rnod_surf[i];
209  }
210 
211  mesh_int_t gmax = global_max(rnod_surf, comm);
212  surfmesh.g_numpts = gmax+1;
213 
214  surfmesh.xyz.resize(surfmesh.l_numpts*3);
215 
216  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
217  T tag = surfmesh.tag[eidx];
218 
219  for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
220  {
221  T g = l2g_surf[surfmesh.con[n]];
222  T old_g = dof2vertex[g];
223 
224  const auto local_idx = g2l_surf[g];
225  const auto orig_idx = g2l_orig[old_g];
226 
227  surfmesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
228  surfmesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
229  surfmesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
230 
231  if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) ||
232  !std::isfinite(surfmesh.xyz[local_idx*3+1]) ||
233  !std::isfinite(surfmesh.xyz[local_idx*3+2]))
234  {
235  fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%lld \n",
236  static_cast<long long>(tag));
237  exit(1);
238  }
239  }
240  }
241 }
242 
243 }
244 #endif
245 #endif
mesh_int_t
opencarp::local_index_t mesh_int_t
Definition: SF_container.h:31
SF_mesh_io.h
Functions related to mesh IO.
backend.h
Starts and stops the linear algebra backend the flavor parameter selects.
EXIT
#define EXIT(A)
finalize the backend, if it runs, and exit
Definition: backend.h:23
SF::vector::size
size_t size() const
The current size of the vector.
Definition: SF_vector.h:89
SF::vector::resize
void resize(size_t n)
Resize a vector.
Definition: SF_vector.h:194
hashmap::unordered_map
Definition: hashmap.hpp:242
hashmap::unordered_map::reserve
void reserve(size_t n)
Definition: hashmap.hpp:719
hashmap::unordered_set< int >
hashmap::unordered_set::find
iterator find(const K &key)
Definition: hashmap.hpp:1081
hashmap::unordered_set::end
iterator end()
Definition: hashmap.hpp:1146
SF
Definition: dense_mat.hpp:19
SF::global_max
T global_max(const vector< T > &vec, MPI_Comm comm)
Compute the global maximum of a distributed vector.
Definition: SF_network.h:141
SF::SF_nbr
SF_nbr
Enumeration encoding the different supported numberings.
Definition: SF_container.h:185
Generated by doxygen 1.9.1