@@ -207,15 +207,14 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
207207 {
208208 THROW_ALGORITHM_PROCESSING_ERROR (" The scalp surface (second module input, field, SCALP_TRI_SURF_MESH) must be a triangle mesh. " );
209209 }
210-
210+
211211 VMesh* mesh_vmesh = mesh->vmesh ();
212212 VMesh::size_type mesh_num_nodes = mesh_vmesh->num_nodes ();
213213 mesh_vmesh->synchronize (Mesh::NODE_LOCATE_E);
214214 DenseMatrixHandle lhs_knows, elc_elem, elc_elem_typ, elc_elem_def, elc_con_imp;
215215 std::vector<double > electrode_sponge_areas;
216216 index_type refnode_number = get (Parameters::refnode).toInt ();
217217 bool getsurf = get (Parameters::GetContactSurface).toBool ();
218-
219218 double normal_dot_product_bound_ = get (Parameters::normal_dot_product_bound).toDouble ();
220219 double identical_node_location_differce = get (Parameters::pointdistancebound).toDouble ();
221220
@@ -280,7 +279,7 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
280279 algo.set (SplitFieldByConnectedRegionAlgo::SortDomainBySize (), false );
281280 algo.set (SplitFieldByConnectedRegionAlgo::SortAscending (), false );
282281 std::vector<FieldHandle> result = algo.run (elc_tri_surf);
283-
282+
284283 if (result.size ()<=0 )
285284 {
286285 THROW_ALGORITHM_PROCESSING_ERROR (" Splitting input mesh into connected regions failed. " );
@@ -311,14 +310,15 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
311310 VMesh* elc_sponge_surf_vmesh = elc_sponge_surf->vmesh ();
312311 VField* elc_sponge_surf_vfld = elc_sponge_surf->vfield ();
313312 std::vector<double > field_values, impedances;
314-
313+
315314 // / map the electrode sponge center (CreateElectrodeCoil) to generated tDCS electrode geometry (Cleaver),
316315 // / this mapping is meant to map the GUI inputs to actual electrode geometry by having a lookup table
317316 int electrode_sponges=elc_sponge_location->nrows ();
318317 DenseMatrixHandle lookup (new DenseMatrix (electrode_sponges, 1 ));
319318 DenseMatrixHandle distances (new DenseMatrix (electrode_sponges, 1 ));
320319 VMesh::Node::index_type didx;
321320 double distance=0 ;
321+
322322 for (long i=0 ;i<electrode_sponges;i++)
323323 {
324324 Point elc ((*elc_sponge_location)(i,0 ),(*elc_sponge_location)(i,1 ),(*elc_sponge_location)(i,2 )),r;
@@ -330,6 +330,7 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
330330 tmp_mesh->synchronize (Mesh::NODE_LOCATE_E);
331331
332332 tmp_mesh->find_closest_node (distance,r,didx,elc);
333+
333334 if (distance<min_dis)
334335 {
335336 min_dis=distance;
@@ -339,7 +340,7 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
339340 (*lookup)(i,0 )=found_index;
340341 (*distances)(i,0 )=min_dis;
341342 }
342-
343+
343344 // / throw error if the geometry is too far away from predicted location - further away than sponge thickness
344345 for (long i=0 ;i<elc_sponge_location->nrows ();i++)
345346 {
@@ -365,7 +366,6 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
365366 DenseMatrixHandle sponge_center_pojected_onto_scalp (new DenseMatrix (result.size (), 3 ));
366367
367368 mesh_scalp_tri_surf->synchronize (Mesh::NODE_LOCATE_E);
368-
369369 for (long i=0 ;i<result.size ();i++)
370370 {
371371 VMesh* tmp_mesh = result[i]->vmesh ();
@@ -408,7 +408,6 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
408408 mesh_scalp_tri_surf->find_closest_node (distance,p,node_ind,q);
409409 (*sponge_center_pojected_onto_scalp_index)(i,0 )=node_ind;
410410 }
411-
412411 // / determine normal of scalp/electrode sponge inferface center
413412 mesh_scalp_tri_surf->synchronize (Mesh::NORMALS_E);
414413 VMesh::Elem::array_type ca;
@@ -432,11 +431,9 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
432431 mesh_elc_tri_surf->synchronize (Mesh::NODE_LOCATE_E);
433432
434433 int nr_elc_sponge_triangles=0 ;
435-
436434 FieldInformation fi3 (" PointCloudMesh" ,0 ," double" );
437435 FieldHandle estimated_sponge_top_center_points=CreateField (fi3);
438436 VMesh* estimated_sponge_top_center_points_vmesh = estimated_sponge_top_center_points->vmesh ();
439-
440437 // / retrieve electrode thickness from last input by using the lookup table
441438 for (int k=0 ;k<sponge_center_pojected_onto_scalp_normal->nrows ();k++)
442439 {
@@ -545,7 +542,6 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
545542
546543 double dot_sp_o1 = Dot (sp_o1, current_scalp_normal_at_elec);
547544 double dot_sp_o2 = Dot (sp_o2, current_scalp_normal_at_elec);
548-
549545 double x=0 ,y=0 ,z=0 ;
550546 if ( (dot_sp_o1<0 && dot_sp_o2<0 ) || (dot_sp_o1>0 && dot_sp_o2>0 ) )
551547 {
@@ -716,7 +712,7 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
716712 THROW_ALGORITHM_PROCESSING_ERROR (" Internal error: internal field definition is " );
717713 }
718714
719- double area=0.0 ; double prev_elc=0 ;
715+ double area=0.0 , prev_elc= 0 ,tmp_fld_val =0 ;
720716 for (VMesh::Elem::index_type l=0 ; l<elc_sponge_surf_vmesh->num_elems (); l++)
721717 {
722718 (*elc_elem_typ)(l,0 )=2 ; // define triangles to incject currents in TDCS simulations
@@ -753,15 +749,14 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
753749 double x3=pos.x (), y3=pos.y (), z3=pos.z ();
754750
755751 (*elc_elem_def)(l,3 )=0 ;
756- (*elc_elem)(l,0 )=field_values[l];
757- (*elc_con_imp)(l,0 )=impedances[l];
752+ (*elc_elem)(l,0 )=field_values[l];
758753
759754 // /compute surface area of electrode/scalp interface
760755 double area_tmp =y1*z2+z1*y3+y2*z3-z2*y3-z1*y2-y1*z3;
761756 double area_tmp1=z1*x2+x1*z3+z2*x3-x2*z3-x1*z2-z1*x3;
762757 double area_tmp2=x1*y2+y1*x3+x2*y3-y2*x3-y1*x2-x1*y3;
763758 double triangle_area=0.5 * sqrt (area_tmp*area_tmp+area_tmp1*area_tmp1+area_tmp2*area_tmp2);
764- double tmp_fld_val=0 ;
759+ tmp_fld_val=0 ;
765760 elc_sponge_surf_vfld->get_value (tmp_fld_val,l);
766761
767762 if (prev_elc!=tmp_fld_val)
@@ -775,8 +770,21 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
775770 area+=triangle_area;
776771 }
777772 }
773+
778774 area*=1e-6 ;
775+
779776 electrode_sponge_areas.push_back (area);
777+
778+ for (VMesh::Elem::index_type l=0 ; l<elc_sponge_surf_vmesh->num_elems (); l++) // / the impedance (Ohm * m^2) that was provided by the user needs to be defined related to the electrode area. For convenience we do that for the user here.
779+ {
780+ elc_sponge_surf_vfld->get_value (tmp_fld_val,l);
781+ if (tmp_fld_val>=electrode_sponge_areas.size ())
782+ {
783+ THROW_ALGORITHM_PROCESSING_ERROR (" Internal ERROR (should not happen): could not scale contact impedance by electrode surface area (index out of bound). " );
784+ }
785+ (*elc_con_imp)(l,0 )=impedances[tmp_fld_val]*electrode_sponge_areas[tmp_fld_val];
786+ }
787+
780788 DenseMatrixHandle selectmatrixind (new DenseMatrix (mesh_vmesh->num_nodes (), 1 )); // /create indeces for SelectSubMatrix
781789 for (long i=0 ;i<mesh_vmesh->num_nodes ();i++)
782790 {
@@ -802,7 +810,7 @@ DenseMatrixHandle SetupTDCSAlgorithm::create_rhs(FieldHandle mesh, FieldHandle e
802810 double temp = elcs_wanted[i].toDouble ();
803811 min_current += temp;
804812 }
805-
813+
806814 if (std::fabs (min_current) > electode_current_summation_bound)
807815 THROW_ALGORITHM_INPUT_ERROR (" Summed electrode current intensities are greater than 1e-6 mA. The sum should be close to 0 mA !!! " );
808816
@@ -815,7 +823,11 @@ DenseMatrixHandle SetupTDCSAlgorithm::create_rhs(FieldHandle mesh, FieldHandle e
815823 VMesh* mesh_elc_tri_surf = elc_tri_surf->vmesh ();
816824 mesh_elc_tri_surf->synchronize (Mesh::NODE_LOCATE_E);
817825 vmesh->synchronize (Mesh::NODE_LOCATE_E);
818-
826+ for (VMesh::Node::index_type l=0 ; l<vmesh->num_nodes (); l++)
827+ {
828+ (*output)(l,0 )=0 ;
829+ }
830+ double min_dis = get (Parameters::pointdistancebound).toDouble ();
819831 for (VMesh::Node::index_type l=0 ; l<mesh_elc_tri_surf->num_nodes (); l++)
820832 {
821833 Point p,q;
@@ -824,8 +836,14 @@ DenseMatrixHandle SetupTDCSAlgorithm::create_rhs(FieldHandle mesh, FieldHandle e
824836 VMesh::Node::index_type ind;
825837 vmesh->find_closest_node (distance,q,ind,p);
826838 (*output)(ind,0 )=elcs_wanted[l].toDouble ()/1000.0 ;
839+ if (min_dis<distance)
840+ {
841+ std::ostringstream ostr4;
842+ ostr4 << " The electrode locations (4th module input) are further away from them mesh as provided by the GUI (min. distance bound). " << std::endl;
843+ remark (ostr4.str ());
844+ }
827845 }
828-
846+
829847 return output;
830848 } else
831849 {
@@ -859,14 +877,12 @@ boost::tuple<DenseMatrixHandle, DenseMatrixHandle, DenseMatrixHandle, DenseMatri
859877{
860878 if (num_of_elc > max_number_of_electrodes) { THROW_ALGORITHM_INPUT_ERROR (" Number of electrodes (>512) given exceeds what is possible " );} // / number of possible electrodes is currently bound to 512 electrodes in default setting
861879 if (num_of_elc < 0 ) { THROW_ALGORITHM_INPUT_ERROR (" Negative number of electrodes given " );}
862-
863880 if (!mesh) THROW_ALGORITHM_INPUT_ERROR (" Input field (mesh) was not allocated " );
864881 if (!scalp_tri_surf) THROW_ALGORITHM_INPUT_ERROR (" Input field (scalp triangle surface) was not allocated " );
865882 if (!elc_tri_surf) THROW_ALGORITHM_INPUT_ERROR (" Input field (electrode triangle surface) was not allocated " );
866883 if (!elc_sponge_location) THROW_ALGORITHM_INPUT_ERROR (" Input field (electrode triangle surface) was not allocated " );
867-
868884 DenseMatrixHandle rhs=create_rhs (mesh, elc_tri_surf, elcs, num_of_elc); // / get the right-hand-side
869-
885+
870886 DenseMatrixHandle lhs_knowns, elc_element, elc_element_typ, elc_element_def, elc_contact_imp, selectmatrixind;
871887 FieldHandle elec_sponge_surf;
872888 std::vector<double > electrode_sponge_areas;
0 commit comments