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 "SF_mesh_io.h"
17 #include "petsc_utils.h"
18 
19 namespace SF {
20 
36 template<class T, class S>
37 inline void insert_points_ptsData_to_dof( const meshdata<T, S> & mesh_original, meshdata<T, S> & mesh,
38  const SF_nbr numbering,
39  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
40  hashmap::unordered_map<T,T> & vertex2ptsdata,
41  hashmap::unordered_map<T,T> & dof2ptsData)
42 {
43  MPI_Comm comm = mesh_original.comm;
44 
45  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
46  hashmap::unordered_map<T,T> g2l_orig;
47  g2l_orig.reserve(rnod_original.size());
48  for(size_t i=0; i<rnod_original.size(); i++){
49  g2l_orig[rnod_original[i]] = i;
50  }
51 
52  const SF::vector<mesh_int_t> & rnod = mesh.get_numbering(numbering);
53  hashmap::unordered_map<T,T> g2l;
54  hashmap::unordered_map<T,T> l2g;
55  g2l.reserve(rnod.size());
56  l2g.reserve(rnod.size());
57  for(size_t i=0; i<rnod.size(); i++){
58  g2l[rnod[i]] = i;
59  l2g[i] = rnod[i];
60  }
61 
62  mesh_int_t gmax = global_max(rnod, comm);
63  mesh.g_numpts = gmax+1;
64 
65  mesh.xyz.resize(mesh.l_numpts*3);
66 
67  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
68  T tag = mesh.tag[eidx];
69 
70  for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
71  {
72  T g = l2g[mesh.con[n]];
73  T old_g = dof2vertex[g];
74 
75  const auto local_idx = g2l[g];
76  const auto orig_idx = g2l_orig[old_g];
77 
78  mesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
79  mesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
80  mesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
81  dof2ptsData[g] = vertex2ptsdata[old_g];
82 
83  if (!std::isfinite(mesh.xyz[local_idx*3+0]) ||
84  !std::isfinite(mesh.xyz[local_idx*3+1]) ||
85  !std::isfinite(mesh.xyz[local_idx*3+2]))
86  {
87  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",
88  static_cast<long long>(tag));
89  EXIT(EXIT_FAILURE);
90  }
91  }
92  }
93 }
94 
110 template<class T, class S>
111 inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh,
112  const SF_nbr numbering,
113  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
114  hashmap::unordered_set<int> & extra_tags,
115  SF::vector<mesh_int_t> & elemTag_surface)
116 {
117  MPI_Comm comm = mesh_original.comm;
118 
119  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
120  hashmap::unordered_map<T,T> g2l_orig;
121  g2l_orig.reserve(rnod_original.size());
122  for(size_t i=0; i<rnod_original.size(); i++){
123  g2l_orig[rnod_original[i]] = i;
124  }
125 
126  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
127 
128  hashmap::unordered_map<T,T> g2l_surf;
129  hashmap::unordered_map<T,T> l2g_surf;
130  g2l_surf.reserve(rnod_surf.size());
131  l2g_surf.reserve(rnod_surf.size());
132  for(size_t i=0; i<rnod_surf.size(); i++){
133  g2l_surf[rnod_surf[i]] = i;
134  l2g_surf[i] = rnod_surf[i];
135  }
136 
137  mesh_int_t gmax = global_max(rnod_surf, comm);
138  surfmesh.g_numpts = gmax+1;
139 
140  surfmesh.xyz.resize(surfmesh.l_numpts*3);
141 
142  elemTag_surface.resize(surfmesh.l_numelem);
143  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
144  T tag = surfmesh.tag[eidx];
145  elemTag_surface[eidx] = 1; // assigned 1 for for extracellular located on extracellular side
146  if(extra_tags.find(tag) == extra_tags.end())
147  elemTag_surface[eidx] = 2; // assigned 2 for for intracellular located on myocyte
148 
149  for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
150  {
151  T g = l2g_surf[surfmesh.con[n]];
152  T old_g = dof2vertex[g];
153 
154  const auto local_idx = g2l_surf[g];
155  const auto orig_idx = g2l_orig[old_g];
156 
157  surfmesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
158  surfmesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
159  surfmesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
160 
161  if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) ||
162  !std::isfinite(surfmesh.xyz[local_idx*3+1]) ||
163  !std::isfinite(surfmesh.xyz[local_idx*3+2]))
164  {
165  fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%lld \n",
166  static_cast<long long>(tag));
167  exit(1);
168  }
169  }
170  }
171 }
172 
184 template<class T, class S>
185 inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh,
186  const SF_nbr numbering,
187  hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex)
188 {
189  MPI_Comm comm = mesh_original.comm;
190 
191  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
192  hashmap::unordered_map<T,T> g2l_orig;
193  g2l_orig.reserve(rnod_original.size());
194  for(size_t i=0; i<rnod_original.size(); i++){
195  g2l_orig[rnod_original[i]] = i;
196  }
197 
198  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
199 
200  hashmap::unordered_map<T,T> g2l_surf;
201  hashmap::unordered_map<T,T> l2g_surf;
202  g2l_surf.reserve(rnod_surf.size());
203  l2g_surf.reserve(rnod_surf.size());
204  for(size_t i=0; i<rnod_surf.size(); i++){
205  g2l_surf[rnod_surf[i]] = i;
206  l2g_surf[i] = rnod_surf[i];
207  }
208 
209  mesh_int_t gmax = global_max(rnod_surf, comm);
210  surfmesh.g_numpts = gmax+1;
211 
212  surfmesh.xyz.resize(surfmesh.l_numpts*3);
213 
214  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
215  T tag = surfmesh.tag[eidx];
216 
217  for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
218  {
219  T g = l2g_surf[surfmesh.con[n]];
220  T old_g = dof2vertex[g];
221 
222  const auto local_idx = g2l_surf[g];
223  const auto orig_idx = g2l_orig[old_g];
224 
225  surfmesh.xyz[local_idx*3+0] = mesh_original.xyz[orig_idx*3+0];
226  surfmesh.xyz[local_idx*3+1] = mesh_original.xyz[orig_idx*3+1];
227  surfmesh.xyz[local_idx*3+2] = mesh_original.xyz[orig_idx*3+2];
228 
229  if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) ||
230  !std::isfinite(surfmesh.xyz[local_idx*3+1]) ||
231  !std::isfinite(surfmesh.xyz[local_idx*3+2]))
232  {
233  fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%lld \n",
234  static_cast<long long>(tag));
235  exit(1);
236  }
237  }
238  }
239 }
240 
241 }
242 #endif
243 #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.
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