Changes in src/Actions/AnalysisAction/PairCorrelationAction.cpp [112b09:623e89]

File:

• src/Actions/AnalysisAction/PairCorrelationAction.cpp (modified) (4 diffs)

src/Actions/AnalysisAction/PairCorrelationAction.cpp

@@ -6 +6 @@
  */
 
-#include "Helpers/MemDebug.hpp"
-
 #include "Actions/AnalysisAction/PairCorrelationAction.hpp"
-#include "CommandLineParser.hpp"
 #include "analysis_correlation.hpp"
+#include "boundary.hpp"
+#include "linkedcell.hpp"
 #include "log.hpp"
 #include "element.hpp"
+#include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "vector.hpp"
 #include "World.hpp"
 
@@ -37 +38 @@
   Dialog *dialog = UIFactory::getInstance().makeDialog();
   int ranges[3] = {1, 1, 1};
+  double BinEnd = 0.;
   double BinStart = 0.;
-  double BinEnd = 0.;
+  double BinWidth = 0.;
+  molecule *Boundary = NULL;
   string outputname;
   string binoutputname;
@@ -44 +47 @@
   ofstream output;
   ofstream binoutput;
-  const element *elemental1;
-  const element *elemental2;
+  std::vector< element *> elements;
+  string type;
+  Vector Point;
+  BinPairMap *binmap = NULL;
+  MoleculeListClass *molecules = World::getInstance().getMolecules();
 
-  dialog->queryElement("elements", &elemental1, MapOfActions::getInstance().getDescription("elements"));
-  dialog->queryElement("elements", &elemental2, MapOfActions::getInstance().getDescription("elements"));
+  // first dialog: Obtain which type of correlation
+  dialog->queryString(NAME, &type, MapOfActions::getInstance().getDescription(NAME));
+  if(dialog->display()) {
+    delete dialog;
+  } else {
+    delete dialog;
+    return Action::failure;
+  }
+
+  // second dialog: Obtain parameters specific to this type
+  dialog = UIFactory::getInstance().makeDialog();
+  if (type == "P")
+    dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
+  if (type == "S")
+    dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));
+  dialog->queryElement("elements", &elements, MapOfActions::getInstance().getDescription("elements"));
   dialog->queryDouble("bin-start", &BinStart, MapOfActions::getInstance().getDescription("bin-start"));
+  dialog->queryDouble("bin-width", &BinWidth, MapOfActions::getInstance().getDescription("bin-width"));
   dialog->queryDouble("bin-end", &BinEnd, MapOfActions::getInstance().getDescription("bin-end"));
   dialog->queryString("output-file", &outputname, MapOfActions::getInstance().getDescription("output-file"));
@@ -59 +80 @@
     binoutput.open(binoutputname.c_str());
     PairCorrelationMap *correlationmap = NULL;
-    if (periodic)
-      correlationmap = PeriodicPairCorrelation(World::getInstance().getMolecules(), elemental1, elemental2, ranges);
-    else
-      correlationmap = PairCorrelation(World::getInstance().getMolecules(), elemental1, elemental2);
-    //OutputCorrelationToSurface(&output, correlationmap);
-    BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
+    if (type == "E") {
+      PairCorrelationMap *correlationmap = NULL;
+      if (periodic)
+        correlationmap = PeriodicPairCorrelation(World::getInstance().getMolecules(), elements, ranges);
+      else
+        correlationmap = PairCorrelation(World::getInstance().getMolecules(), elements);
+      //OutputCorrelationToSurface(&output, correlationmap);
+      binmap = BinData( correlationmap, BinWidth, BinStart, BinEnd );
+    } else if (type == "P")  {
+      cout << "Point to correlate to is  " << Point << endl;
+      CorrelationToPointMap *correlationmap = NULL;
+      if (periodic)
+        correlationmap  = PeriodicCorrelationToPoint(molecules, elements, &Point, ranges);
+      else
+        correlationmap = CorrelationToPoint(molecules, elements, &Point);
+      //OutputCorrelationToSurface(&output, correlationmap);
+      binmap = BinData( correlationmap, BinWidth, BinStart, BinEnd );
+    } else if (type == "S") {
+      ASSERT(Boundary != NULL, "No molecule specified for SurfaceCorrelation.");
+      const double radius = 4.;
+      double LCWidth = 20.;
+      if (BinEnd > 0) {
+        if (BinEnd > 2.*radius)
+            LCWidth = BinEnd;
+        else
+          LCWidth = 2.*radius;
+      }
+
+      // get the boundary
+      class Tesselation *TesselStruct = NULL;
+      const LinkedCell *LCList = NULL;
+      // find biggest molecule
+      int counter  = molecules->ListOfMolecules.size();
+      bool *Actives = new bool[counter];
+      counter = 0;
+      for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
+        Actives[counter++] = (*BigFinder)->ActiveFlag;
+        (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
+      }
+      LCList = new LinkedCell(Boundary, LCWidth);
+      FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
+      CorrelationToSurfaceMap *surfacemap = NULL;
+      if (periodic)
+        surfacemap = PeriodicCorrelationToSurface( molecules, elements, TesselStruct, LCList, ranges);
+      else
+        surfacemap = CorrelationToSurface( molecules, elements, TesselStruct, LCList);
+      OutputCorrelationToSurface(&output, surfacemap);
+      // check whether radius was appropriate
+      {
+        double start; double end;
+        GetMinMax( surfacemap, start, end);
+        if (LCWidth < end)
+          DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl);
+      }
+      binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
+    } else
+      return Action::failure;
     OutputCorrelation ( &binoutput, binmap );
     output.close();