Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FilterElectrogram.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterElectrogram.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Outputs virtual electrograms at specific locations.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <iomanip>
39 
40 namespace Nektar
41 {
42  std::string FilterElectrogram::className =
44  "Electrogram",
46 
47  /**
48  *
49  */
52  const std::map<std::string, std::string> &pParams)
53  : Filter(pSession)
54  {
55  if (pParams.find("OutputFile") == pParams.end())
56  {
57  m_outputFile = m_session->GetSessionName();
58  }
59  else
60  {
61  ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
62  "Missing parameter 'OutputFile'.");
63  m_outputFile = pParams.find("OutputFile")->second;
64  }
65  if (!(m_outputFile.length() >= 4
66  && m_outputFile.substr(m_outputFile.length() - 4) == ".ecg"))
67  {
68  m_outputFile += ".ecg";
69  }
70 
71  if (pParams.find("OutputFrequency") == pParams.end())
72  {
74  }
75  else
76  {
78  atoi(pParams.find("OutputFrequency")->second.c_str());
79  }
80 
81  ASSERTL0(pParams.find("Points") != pParams.end(),
82  "Missing parameter 'Points'.");
83  m_electrogramStream.str(pParams.find("Points")->second);
84  m_index = 0;
85  }
86 
87 
88  /**
89  *
90  */
92  {
93 
94  }
95 
96 
97  /**
98  *
99  */
101  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
102  const NekDouble &time)
103  {
105  "No history points in stream.");
106 
107  m_index = 0;
108  Array<OneD, NekDouble> gloCoord(3,0.0);
109  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
110 
111  // Read electrogram points
112  // Always use dim = 3 to allow electrode to be above surface
113  const int dim = 3;
114  int i = 0;
115 
116  while (!m_electrogramStream.fail())
117  {
118  m_electrogramStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
119  if (!m_electrogramStream.fail())
120  {
123  ::AllocateSharedPtr(dim, i, gloCoord[0],
124  gloCoord[1], gloCoord[2]);
125 
126  m_electrogramPoints.push_back(vert);
127  ++i;
128  }
129  }
130 
131  if (vComm->GetRank() == 0)
132  {
133  // Open output stream
134  m_outputStream.open(m_outputFile.c_str());
135  m_outputStream << "# Electrogram data for variables (:";
136 
137  for (i = 0; i < pFields.num_elements(); ++i)
138  {
139  m_outputStream << m_session->GetVariable(i) <<",";
140  }
141 
142  m_outputStream << ") at points:" << endl;
143 
144  for (i = 0; i < m_electrogramPoints.size(); ++i)
145  {
146  m_electrogramPoints[i]->GetCoords( gloCoord[0],
147  gloCoord[1],
148  gloCoord[2]);
149 
150  m_outputStream << "# \t" << i;
151  m_outputStream.width(8);
152  m_outputStream << gloCoord[0];
153  m_outputStream.width(8);
154  m_outputStream << gloCoord[1];
155  m_outputStream.width(8);
156  m_outputStream << gloCoord[2];
157  m_outputStream << endl;
158  }
159  }
160 
161  // Compute the distance function for each electrogram point
162  const unsigned int nq = pFields[0]->GetNpoints();
163  const unsigned int npts = m_electrogramPoints.size();
164  NekDouble px, py, pz;
165  m_grad_R_x = Array<OneD, Array<OneD, NekDouble> >(npts);
166  m_grad_R_y = Array<OneD, Array<OneD, NekDouble> >(npts);
167  m_grad_R_z = Array<OneD, Array<OneD, NekDouble> >(npts);
168 
169  Array<OneD, NekDouble> x(nq);
170  Array<OneD, NekDouble> y(nq);
171  Array<OneD, NekDouble> z(nq);
172 
173  Array<OneD, NekDouble> oneOverR(nq);
174  for (unsigned int i = 0; i < npts; ++i)
175  {
176  m_grad_R_x[i] = Array<OneD, NekDouble>(nq);
177  m_grad_R_y[i] = Array<OneD, NekDouble>(nq);
178  m_grad_R_z[i] = Array<OneD, NekDouble>(nq);
179 
180  // Compute 1/R
181  m_electrogramPoints[i]->GetCoords(px,py,pz);
182 
183  pFields[0]->GetCoords(x,y,z);
184 
185  Vmath::Sadd (nq, -px, x, 1, x, 1);
186  Vmath::Sadd (nq, -py, y, 1, y, 1);
187  Vmath::Sadd (nq, -pz, z, 1, z, 1);
188  Vmath::Vvtvvtp(nq, x, 1, x, 1, y, 1, y, 1, oneOverR, 1);
189  Vmath::Vvtvp (nq, z, 1, z, 1, oneOverR, 1, oneOverR, 1);
190  Vmath::Vsqrt (nq, oneOverR, 1, oneOverR, 1);
191  Vmath::Sdiv (nq, 1.0, oneOverR, 1, oneOverR, 1);
192 
193  // Compute grad 1/R
194  pFields[0]->PhysDeriv(oneOverR, m_grad_R_x[i], m_grad_R_y[i],
195  m_grad_R_z[i]);
196  }
197 
198  // Compute electrogram point for initial condition
199  v_Update(pFields, time);
200  }
201 
202 
203  /**
204  *
205  */
207  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
208  const NekDouble &time)
209  {
210  // Only output every m_outputFrequency.
211  if ((m_index++) % m_outputFrequency)
212  {
213  return;
214  }
215 
216  const unsigned int nq = pFields[0]->GetNpoints();
217  const unsigned int npoints = m_electrogramPoints.size();
218  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
219 
220  unsigned int i = 0;
221  Array<OneD, NekDouble> e(npoints);
222 
223  // Compute grad V
224  Array<OneD, NekDouble> grad_V_x(nq), grad_V_y(nq), grad_V_z(nq);
225  pFields[0]->PhysDeriv(pFields[0]->GetPhys(),
226  grad_V_x, grad_V_y, grad_V_z);
227 
228  for (i = 0; i < npoints; ++i)
229  {
230  // Multiply together
231  Array<OneD, NekDouble> output(nq);
232  Vmath::Vvtvvtp(nq, m_grad_R_x[i], 1, grad_V_x, 1, m_grad_R_y[i], 1,
233  grad_V_y, 1, output, 1);
234  Vmath::Vvtvp (nq, m_grad_R_z[i], 1, grad_V_z, 1, output, 1,
235  output, 1);
236 
237  e[i] = pFields[0]->Integral(output);
238  }
239 
240  // Exchange history data
241  // This could be improved to reduce communication but works for now
242  vComm->AllReduce(e, LibUtilities::ReduceSum);
243 
244  // Only the root process writes out electrogram data
245  if (vComm->GetRank() == 0)
246  {
247  m_outputStream.width(8);
248  m_outputStream << setprecision(6) << time;
249 
250  // Write data values point by point
251  for (i = 0; i < m_electrogramPoints.size(); ++i)
252  {
253  m_outputStream.width(25);
254  m_outputStream << setprecision(16) << e[i];
255  }
256  m_outputStream << endl;
257  }
258  }
259 
260 
261  /**
262  *
263  */
265  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
266  const NekDouble &time)
267  {
268  if (pFields[0]->GetComm()->GetRank() == 0)
269  {
270  m_outputStream.close();
271  }
272  }
273 
274 
275  /**
276  *
277  */
279  {
280  return true;
281  }
282 }