Nektar++
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 ParamMap &pParams)
53  : Filter(pSession)
54 {
55  ParamMap::const_iterator it;
56 
57  // OutputFile
58  it = pParams.find("OutputFile");
59  if (it == pParams.end())
60  {
61  m_outputFile = m_session->GetSessionName();
62  }
63  else
64  {
65  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
66  m_outputFile = it->second;
67  }
68  if (!(m_outputFile.length() >= 4
69  && m_outputFile.substr(m_outputFile.length() - 4) == ".ecg"))
70  {
71  m_outputFile += ".ecg";
72  }
73 
74  // OutputFrequency
75  it = pParams.find("OutputFrequency");
76  if (it == pParams.end())
77  {
79  }
80  else
81  {
82  LibUtilities::Equation equ(m_session, it->second);
83  m_outputFrequency = floor(equ.Evaluate());
84  }
85 
86  // Points
87  it = pParams.find("Points");
88  ASSERTL0(it != pParams.end(), "Missing parameter 'Points'.");
89  m_electrogramStream.str(it->second);
90  m_index = 0;
91 }
92 
93 
94 /**
95  *
96  */
98 {
99 
100 }
101 
102 
103 /**
104  *
105  */
108  const NekDouble &time)
109 {
111  "No history points in stream.");
112 
113  m_index = 0;
114  Array<OneD, NekDouble> gloCoord(3,0.0);
115  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
116 
117  // Read electrogram points
118  // Always use dim = 3 to allow electrode to be above surface
119  const int dim = 3;
120  int i = 0;
121 
122  while (!m_electrogramStream.fail())
123  {
124  m_electrogramStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
125  if (!m_electrogramStream.fail())
126  {
129  ::AllocateSharedPtr(dim, i, gloCoord[0],
130  gloCoord[1], gloCoord[2]);
131 
132  m_electrogramPoints.push_back(vert);
133  ++i;
134  }
135  }
136 
137  if (vComm->GetRank() == 0)
138  {
139  // Open output stream
140  m_outputStream.open(m_outputFile.c_str());
141  m_outputStream << "# Electrogram data for variables (:";
142 
143  for (i = 0; i < pFields.num_elements(); ++i)
144  {
145  m_outputStream << m_session->GetVariable(i) <<",";
146  }
147 
148  m_outputStream << ") at points:" << endl;
149 
150  for (i = 0; i < m_electrogramPoints.size(); ++i)
151  {
152  m_electrogramPoints[i]->GetCoords( gloCoord[0],
153  gloCoord[1],
154  gloCoord[2]);
155 
156  m_outputStream << "# \t" << i;
157  m_outputStream.width(8);
158  m_outputStream << gloCoord[0];
159  m_outputStream.width(8);
160  m_outputStream << gloCoord[1];
161  m_outputStream.width(8);
162  m_outputStream << gloCoord[2];
163  m_outputStream << endl;
164  }
165  }
166 
167  // Compute the distance function for each electrogram point
168  const unsigned int nq = pFields[0]->GetNpoints();
169  const unsigned int npts = m_electrogramPoints.size();
170  NekDouble px, py, pz;
174 
178 
179  Array<OneD, NekDouble> oneOverR(nq);
180  for (unsigned int i = 0; i < npts; ++i)
181  {
183  m_grad_R_y[i] = Array<OneD, NekDouble>(nq);
184  m_grad_R_z[i] = Array<OneD, NekDouble>(nq);
185 
186  // Compute 1/R
187  m_electrogramPoints[i]->GetCoords(px,py,pz);
188 
189  pFields[0]->GetCoords(x,y,z);
190 
191  Vmath::Sadd (nq, -px, x, 1, x, 1);
192  Vmath::Sadd (nq, -py, y, 1, y, 1);
193  Vmath::Sadd (nq, -pz, z, 1, z, 1);
194  Vmath::Vvtvvtp(nq, x, 1, x, 1, y, 1, y, 1, oneOverR, 1);
195  Vmath::Vvtvp (nq, z, 1, z, 1, oneOverR, 1, oneOverR, 1);
196  Vmath::Vsqrt (nq, oneOverR, 1, oneOverR, 1);
197  Vmath::Sdiv (nq, 1.0, oneOverR, 1, oneOverR, 1);
198 
199  // Compute grad 1/R
200  pFields[0]->PhysDeriv(oneOverR, m_grad_R_x[i], m_grad_R_y[i],
201  m_grad_R_z[i]);
202  }
203 
204  // Compute electrogram point for initial condition
205  v_Update(pFields, time);
206 }
207 
208 
209 /**
210  *
211  */
214  const NekDouble &time)
215 {
216  // Only output every m_outputFrequency.
217  if ((m_index++) % m_outputFrequency)
218  {
219  return;
220  }
221 
222  const unsigned int nq = pFields[0]->GetNpoints();
223  const unsigned int npoints = m_electrogramPoints.size();
224  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
225 
226  unsigned int i = 0;
227  Array<OneD, NekDouble> e(npoints);
228 
229  // Compute grad V
230  Array<OneD, NekDouble> grad_V_x(nq), grad_V_y(nq), grad_V_z(nq);
231  pFields[0]->PhysDeriv(pFields[0]->GetPhys(),
232  grad_V_x, grad_V_y, grad_V_z);
233 
234  for (i = 0; i < npoints; ++i)
235  {
236  // Multiply together
237  Array<OneD, NekDouble> output(nq);
238  Vmath::Vvtvvtp(nq, m_grad_R_x[i], 1, grad_V_x, 1, m_grad_R_y[i], 1,
239  grad_V_y, 1, output, 1);
240  Vmath::Vvtvp (nq, m_grad_R_z[i], 1, grad_V_z, 1, output, 1,
241  output, 1);
242 
243  e[i] = pFields[0]->Integral(output);
244  }
245 
246  // Exchange history data
247  // This could be improved to reduce communication but works for now
248  vComm->AllReduce(e, LibUtilities::ReduceSum);
249 
250  // Only the root process writes out electrogram data
251  if (vComm->GetRank() == 0)
252  {
253  m_outputStream.width(8);
254  m_outputStream << setprecision(6) << time;
255 
256  // Write data values point by point
257  for (i = 0; i < m_electrogramPoints.size(); ++i)
258  {
259  m_outputStream.width(25);
260  m_outputStream << setprecision(16) << e[i];
261  }
262  m_outputStream << endl;
263  }
264 }
265 
266 
267 /**
268  *
269  */
272  const NekDouble &time)
273 {
274  if (pFields[0]->GetComm()->GetRank() == 0)
275  {
276  m_outputStream.close();
277  }
278 }
279 
280 
281 /**
282  *
283  */
285 {
286  return true;
287 }
288 }
unsigned int m_index
Counts number of calls to update (number of timesteps)
std::string m_outputFile
Filename to output electrogram data to.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Array< OneD, Array< OneD, NekDouble > > m_grad_R_y
Gradient of the radius from each electrogram point in y-direction.
std::stringstream m_electrogramStream
Point coordinate input string.
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Finalise the electrogram filter and close output file.
std::ofstream m_outputStream
Output file stream for electrogram data.
virtual bool v_IsTimeDependent()
Filter is time-dependent and should be called at each time-step.
SpatialDomains::PointGeomVector m_electrogramPoints
List of electrogram points.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
virtual void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Compute extracellular potential at egm points at current time.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
static std::string className
Name of the class.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
unsigned int m_outputFrequency
Number of timesteps between outputs.
static SolverUtils::FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
Creates an instance of this class.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
NekDouble Evaluate() const
Definition: Equation.h:102
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
FilterElectrogram(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
Electrogram filter constructor.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:523
Array< OneD, Array< OneD, NekDouble > > m_grad_R_z
Gradient of the radius from each electrogram point in z-direction.
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
virtual ~FilterElectrogram()
Electrogram filter destructor.
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Initialises the electrogram filter and open output file.
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
Array< OneD, Array< OneD, NekDouble > > m_grad_R_x
Gradient of the radius from each electrogram point in x-direction.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215