small update

Author: moritzdannhauer

src/Core/Algorithms/Field/InterfaceWithCleaverAlgorithm.cc

@@ -31,7 +31,10 @@ ToDo: Padding is always enabled because of exit() in cleaver lib
 */
 
 ///TODO: fix include path to remove Externals/ part
+
 #include <Externals/cleaver/lib/FloatField.h>
+#include <Externals/cleaver/lib/vec3.h>
+#include <Externals/cleaver/lib/BoundingBox.h>
 #include <Externals/cleaver/lib/Cleaver.h>
 #include <Externals/cleaver/lib/InverseField.h>
 #include <Externals/cleaver/lib/PaddedVolume.h>
@@ -49,6 +52,7 @@ ToDo: Padding is always enabled because of exit() in cleaver lib
 #include <Core/Utils/StringUtil.h>
 #include <boost/scoped_ptr.hpp>
 #include <Core/Logging/Log.h>
+#include <Core/Math/MiscMath.h>
 
 using namespace SCIRun::Core::Algorithms;
 using namespace SCIRun::Core::Algorithms::Fields;
@@ -74,24 +78,38 @@ InterfaceWithCleaverAlgorithm::InterfaceWithCleaverAlgorithm()
   addParameter(VolumeScalingZ,1.0);  
 }
 
-boost::shared_ptr<Cleaver::ScalarField> InterfaceWithCleaverAlgorithm::makeCleaverFieldFromLatVol(FieldHandle field)
+boost::shared_ptr<Cleaver::ScalarField> InterfaceWithCleaverAlgorithm::makeCleaverFieldFromLatVol(FieldHandle field )
 {
   //TODO: this function assumes input is completely valid, may want to move various checks from run() function here.
-
   VMesh*  vmesh   = field->vmesh();
   VField* vfield = field->vfield();
   VMesh::dimension_type dims;
   vmesh->get_dimensions( dims ); 
+  
   float* ptr = static_cast<float*>(vfield->fdata_pointer());
-  auto cleaverField = boost::make_shared<Cleaver::FloatField>(dims[0], dims[1], dims[2], ptr);
-
-  //0 = constant, 1 = linear
-  if (0 == vfield->basis_order())
-    cleaverField->setCenter(Cleaver::FloatField::CellCentered);
-  else if (1 == vfield->basis_order())
-    cleaverField->setCenter(Cleaver::FloatField::NodeCentered);
-  //else: TODO? handle other bases, probably by throwing exception.
-
+  
+  auto cleaverField = boost::make_shared<Cleaver::FloatField>(dims[0], dims[1], dims[2], ptr);  
+    
+  cleaverField->setCenter(Cleaver::FloatField::CellCentered);
+    
+  BBox bbox=vmesh->get_bounding_box();
+  Point bmin, bmax;
+  if (bbox.valid())
+  {
+   bmin = bbox.min();
+   bmax = bbox.max();
+  }
+  Cleaver::BoundingBox  bb = Cleaver::BoundingBox(Cleaver::vec3::zero, Cleaver::vec3(dims[0],dims[1],dims[2]));
+  cleaverField->setBounds(bb);
+  const Transform &transform = vmesh->get_transform();
+  
+  int x_spacing=fabs(transform.get_mat_val(0,0)), y_spacing=fabs(transform.get_mat_val(1,1)), z_spacing=fabs(transform.get_mat_val(2,2));
+  if (IsNan(x_spacing)) x_spacing=1;
+  if (IsNan(y_spacing)) y_spacing=1;
+  if (IsNan(z_spacing)) z_spacing=1;
+  
+  cleaverField->setScale(Cleaver::vec3(x_spacing,y_spacing,z_spacing));
+  
   //TODO: Cleaver FloatField needs more setup (constructor is not sufficient)
   // 1.
   // setBounds(BoundingBox). Convert vmesh->get_bounding_box() to Cleaver BoundingBox.
@@ -105,7 +123,6 @@ boost::shared_ptr<Cleaver::ScalarField> InterfaceWithCleaverAlgorithm::makeCleav
 FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& input) const
 {
   FieldHandle output;      
-
   std::vector<FieldHandle> inputs;
   std::copy_if(input.begin(), input.end(), std::back_inserter(inputs), [](FieldHandle f) { return f; });
 
@@ -119,10 +136,11 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
     THROW_ALGORITHM_INPUT_ERROR(" At least 2 indicator functions stored as float values are needed to run cleaver! " );
     return FieldHandle();
   }
-
+  
   std::vector<boost::shared_ptr<Cleaver::ScalarField>> fields;  
+
   VMesh::dimension_type dims; int x=0,y=0,z=0; 
-  for (size_t p=1; p<inputs.size(); p++)
+  for (size_t p=0; p<inputs.size(); p++)
   {
     FieldHandle input = inputs[p];
     VMesh*  imesh1   = input->vmesh();
@@ -141,7 +159,7 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
       }
 
       imesh1->get_dimensions( dims ); 
-      if (p==1)
+      if (p==0)
       {
         x=dims[0]; y=dims[1]; z=dims[2];
         if (x<1 || y<1 || z<1)
@@ -163,26 +181,39 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
         return FieldHandle();
       }    
 
+      //0 = constant, 1 = linear
+       if (0 != vfield1->basis_order())
+       {
+        THROW_ALGORITHM_INPUT_ERROR("Input data need to be defined on elements. You can use the module MapFieldDataFromNodeToElem to do that.");
+       } 
+      
       if (vfield1->is_float())
-      {
+      {    
         float* ptr = static_cast<float*>(vfield1->fdata_pointer());
-        if (ptr)
-        {
+	if (ptr)
+        {	
           fields.push_back(makeCleaverFieldFromLatVol(input));
         }
         else
         {
           THROW_ALGORITHM_INPUT_ERROR(" float field is NULL pointer");
           return FieldHandle();
         }
+      } else
+      {
+       THROW_ALGORITHM_INPUT_ERROR(" Input field needs to be a structured mesh (LatVOL) with float values defnied on the elements. "); 
       }
 
     }
 
   }
 
   boost::shared_ptr<Cleaver::Volume> volume(new Cleaver::Volume(toVectorOfRawPointers(fields)));
-
+/*
+  Cleaver::BoundingBox m_bounds=fields[0]->bounds();
+  std::cout << "bsize: " << m_bounds.size.x << " " << m_bounds.size.y << " " <<  m_bounds.size.z << std::endl;
+  std::cout << "borin: " << m_bounds.origin.x << " " << m_bounds.origin.y << " " << m_bounds.origin.z << std::endl;
+  */ 
   const double xScale = get(VolumeScalingX).toDouble();
   const double yScale = get(VolumeScalingY).toDouble();
   const double zScale = get(VolumeScalingZ).toDouble();
@@ -199,28 +230,31 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
       double newX = xScale*volume->size().x;
       double newY = yScale*volume->size().y;
       double newZ = zScale*volume->size().z;
-      std::cout << "Cleaver setting volume size to " << newX << "," << newY << "," << newZ << std::endl;
       volume->setSize(newX, newY, newZ);
     }
     else // None
+    {
       volume->setSize(dims[0],dims[1],dims[2]);
+    }
   }
   else
   {
     THROW_ALGORITHM_INPUT_ERROR(" Invalid Scaling. Use Input sizes.");
   }
-
+    
   /// Padding is now optional! 
   boost::shared_ptr<Cleaver::AbstractVolume> paddedVolume(volume);
-
   const bool verbose = get(Verbose).toBool();
   const bool pad = get(Padding).toBool();
   if (pad)
     paddedVolume.reset(new Cleaver::PaddedVolume(volume.get()));
+    
+  std::cout << "Creating Mesh with Volume Size " << paddedVolume->size().toString() << std::endl;
+    
   boost::scoped_ptr<Cleaver::TetMesh> mesh(Cleaver::createMeshFromVolume(paddedVolume.get(), verbose));
 
   FieldInformation fi("TetVolMesh",0,"double");   ///create output field
-
+  
   output = CreateField(fi);
   auto omesh = output->vmesh();
   auto ofield = output->vfield();
@@ -247,16 +281,15 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
     vdata[2]=mesh->tets[i]->verts[2]->tm_v_index;
     vdata[3]=mesh->tets[i]->verts[3]->tm_v_index;
     omesh->add_elem(vdata);
-    auto mat_label = mesh->tets[i]->mat_label+1;
+    auto mat_label = mesh->tets[i]->mat_label;
     values[i]=mat_label;
   }
   ofield->resize_values();
   ofield->set_values(values);
   mesh->computeAngles();
   std::ostringstream ostr1;
-  ostr1 << "Number of tetrahedral elements:" << ofield->vmesh()->num_elems() <<  std::endl;
-  ostr1 << "Number of tetrahedral nodes:" << ofield->vmesh()->num_nodes() << std::endl;
-  ostr1 << "Number of tetrahedral nodes:" << ofield->vmesh()->num_nodes() << std::endl;
+  ostr1 << "Number of tetrahedral elements:" << nr_of_tets <<  std::endl;
+  ostr1 << "Number of tetrahedral nodes:" << nr_of_verts << std::endl;
   ostr1 << "Worst Angle (min):" <<  mesh->min_angle << std::endl;
   ostr1 << "Worst Angle (max):" <<  mesh->max_angle << std::endl;
   ostr1 << "Volume:" << volume->size().toString() << std::endl;