ModifyComponent.cpp

// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.

#include <vector>
#include <algorithm>
#include "GmshConfig.h"
#include "ModifyComponent.h"
#include "OctreePost.h"
#include "mathEvaluator.h"

StringXNumber ModifyComponentOptions_Number[] = {
  {GMSH_FULLRC, "Component", NULL, -1.},
  {GMSH_FULLRC, "TimeStep", NULL, -1.},
  {GMSH_FULLRC, "View", NULL, -1.},
  {GMSH_FULLRC, "OtherTimeStep", NULL, -1.},
  {GMSH_FULLRC, "OtherView", NULL, -1.},
  {GMSH_FULLRC, "ForceInterpolation", NULL, 0.}
};

StringXString ModifyComponentOptions_String[] = {
  {GMSH_FULLRC, "Expression", NULL, "v0 * Sin(x)"},
  {GMSH_FULLRC, "Expression1", NULL, ""},
  {GMSH_FULLRC, "Expression2", NULL, ""},
  {GMSH_FULLRC, "Expression3", NULL, ""},
  {GMSH_FULLRC, "Expression4", NULL, ""},
  {GMSH_FULLRC, "Expression5", NULL, ""},
  {GMSH_FULLRC, "Expression6", NULL, ""},
  {GMSH_FULLRC, "Expression7", NULL, ""},
  {GMSH_FULLRC, "Expression8", NULL, ""}
};

extern "C"
{
  GMSH_Plugin *GMSH_RegisterModifyComponentPlugin()
  {
    return new GMSH_ModifyComponentPlugin();
  }
}

std::string GMSH_ModifyComponentPlugin::getHelp() const
{
  return "Plugin(ModifyComponent) sets the `Component'-th "
    "component of the `TimeStep'-th time step in the "
    "view `View' to the expression `Expression'.\n\n"
    "`Expression' can contain:\n\n"
    "- the usual mathematical functions (Log, Sqrt, "
    "Sin, Cos, Fabs, ...) and operators (+, -, *, /, ^);\n\n"
    "If only `Expression' is given (and `Expression1', "
    "..., `Expression8' are all empty), the plugin "
    "creates a scalar view. If `Expression', `Expression1' "
    "and/or `Expression2' are given (and `Expression3', "
    "..., `Expression8' are all empty) the plugin creates "
    "a vector view. Otherwise the plugin creates a tensor "
    "view.\n\n"
    "- the symbols x, y and z, to retrieve the "
    "coordinates of the current node;\n\n"
    "- the symbols Time and TimeStep, to retrieve the "
    "current time and time step values;\n\n"
    "- the symbol v, to retrieve the `Component'-th "
    "component of the field in `View' at the "
    "`TimeStep'-th time step;\n\n"
    "- the symbols v0, v1, v2, ..., v8, to retrieve each "
    "component of the field in `View' at the "
    "`TimeStep'-th time step;\n\n"
    "- the symbol w, to retrieve the `Component'-th "
    "component of the field in `OtherView' at the "
    "`OtherTimeStep'-th time step. If `OtherView' "
    "and `View' are based on different spatial grids, "
    "or if their data types are different, `OtherView' "
    "is interpolated onto `View';\n\n"
    "- the symbols w0, w1, w2, ..., w8, to retrieve each "
    "component of the field in `OtherView' at the "
    "`OtherTimeStep'-th time step.\n\n"
    "If `TimeStep' < 0, the plugin automatically loops "
    "over all the time steps in `View' and evaluates "
    "`Expression' for each one.\n\n"
    "If `OtherTimeStep' < 0, the plugin uses `TimeStep' "
    "instead.\n\n"
    "If `Component' < 0, the plugin automatically  ops\n"
    "over all the components in the view and "
    "evaluates `Expression' for each one.\n\n"
    "If `View' < 0, the plugin is run on the current view.\n\n"
    "If `OtherView' < 0, the plugin uses `View' instead.\n\n"
    "Plugin(ModifyComponent) is executed in-place.";
}

int GMSH_ModifyComponentPlugin::getNbOptions() const
{
  return sizeof(ModifyComponentOptions_Number) / sizeof(StringXNumber);
}

StringXNumber *GMSH_ModifyComponentPlugin::getOption(int iopt)
{
  return &ModifyComponentOptions_Number[iopt];
}

int GMSH_ModifyComponentPlugin::getNbOptionsStr() const
{
  return sizeof(ModifyComponentOptions_String) / sizeof(StringXString);
}

StringXString *GMSH_ModifyComponentPlugin::getOptionStr(int iopt)
{
  return &ModifyComponentOptions_String[iopt];
}

PView *GMSH_ModifyComponentPlugin::execute(PView *view)
{
  int component = (int)ModifyComponentOptions_Number[0].def;
  int timeStep = (int)ModifyComponentOptions_Number[1].def;
  int iView = (int)ModifyComponentOptions_Number[2].def;
  int otherTimeStep = (int)ModifyComponentOptions_Number[3].def;
  int otherView = (int)ModifyComponentOptions_Number[4].def;
  int forceInterpolation = (int)ModifyComponentOptions_Number[5].def;

  PView *v1 = getView(iView, view);
  if(!v1) return view;

  PViewData *data1 = v1->getData();

  if(timeStep > data1->getNumTimeSteps() - 1){
    Msg::Error("Invalid time step (%d) in View[%d]: using step 0 instead",
               timeStep, v1->getIndex());
    timeStep = 0;
  }

  PView *v2 = v1;

  if(otherView >= 0){
    if(otherView < (int)PView::list.size())
      v2 = PView::list[otherView];
    else
      Msg::Error("View[%d] does not exist: using self", otherView);
  }

  PViewData *data2 = getPossiblyAdaptiveData(v2);

  if(otherTimeStep < 0 && data2->getNumTimeSteps() != data1->getNumTimeSteps()){
    Msg::Error("Number of time steps don't match: using step 0");
    otherTimeStep = 0;
  }
  else if(otherTimeStep > data2->getNumTimeSteps() - 1){
    Msg::Error("Invalid time step (%d) in View[%d]: using step 0 instead",
               otherTimeStep, v2->getIndex());
    otherTimeStep = 0;
  }

  std::vector<std::string> expressions(9);
  for(int i = 0; i < 9; i++) expressions[i] = ModifyComponentOptions_String[i].def;
  int numComp3;
  if(expressions[3].size() || expressions[4].size() || expressions[5].size() ||
     expressions[6].size() || expressions[7].size() || expressions[8].size()){
    numComp3 = 9;
    for(int i = 0; i < 9; i++)
      if(expressions[i].empty()) expressions[i] = "0";
  }
  else if(expressions[1].size() || expressions[2].size()){
    numComp3 = 3;
    for(int i = 0; i < 3; i++)
      if(expressions[i].empty()) expressions[i] = "0";
  }
  else{
    numComp3 = 1;
  }
  expressions.resize(numComp3);

  const char *names[] =
    {"x", "y", "z", "Time", "TimeStep",
     "v", "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8",
     "w", "w0", "w1", "w2", "w3", "w4", "w5", "w6", "w7", "w8"};
  unsigned int numVariables = sizeof(names) / sizeof(names[0]);
  std::vector<std::string> variables(numVariables);
  for(unsigned int i = 0; i < numVariables; i++) variables[i] = names[i];
  mathEvaluator f(expressions, variables);
  if(expressions.empty()) return view;
  std::vector<double> values(numVariables), res(numComp3);

  OctreePost *octree = 0;
  if(forceInterpolation ||
     (data1->getNumEntities() != data2->getNumEntities()) ||
     (data1->getNumElements() != data2->getNumElements())){
    Msg::Info("Other view based on different grid: interpolating...");
    octree = new OctreePost(v2);
  }

  for(int step = 0; step < data1->getNumTimeSteps(); step++){
    if(timeStep >= 0 && timeStep != step) continue;

    double time = data1->getTime(step);
    int step2 = (otherTimeStep < 0) ? step : otherTimeStep;

    if(data1->isNodeData()){
      // tag all the nodes with "0" (the default tag)
      for(int ent = 0; ent < data1->getNumEntities(step); ent++){
        for(int ele = 0; ele < data1->getNumElements(step, ent); ele++){
          if(data1->skipElement(step, ent, ele)) continue;
          for(int nod = 0; nod < data1->getNumNodes(step, ent, ele); nod++)
            data1->tagNode(step, ent, ele, nod, 0);
        }
      }
    }

    for(int ent = 0; ent < data1->getNumEntities(step); ent++){
      for(int ele = 0; ele < data1->getNumElements(step, ent); ele++){
        if(data1->skipElement(step, ent, ele)) continue;
        int numComp = data1->getNumComponents(step, ent, ele);
        int numComp2 = octree ? 9 : data2->getNumComponents(step2, ent, ele);
        int numNodes = data1->getNumNodes(step, ent, ele);
        std::vector<double> x(numNodes), y(numNodes), z(numNodes);
        std::vector<int> tag(numNodes);
        for(int nod = 0; nod < numNodes; nod++)
          tag[nod] = data1->getNode(step, ent, ele, nod, x[nod], y[nod], z[nod]);
        for(int nod = 0; nod < numNodes; nod++){
          if(data1->isNodeData() && tag[nod]) continue; // node has already been modified
          std::vector<double> v(std::max(9, numComp), 0.);
          for(int comp = 0; comp < numComp; comp++)
            data1->getValue(step, ent, ele, nod, comp, v[comp]);
          std::vector<double> w(std::max(9, numComp2), 0.);
          if(octree){
            int qn = forceInterpolation ? numNodes : 0;
            if(!octree->searchScalar(x[nod], y[nod], z[nod], &w[0], step2,
                                     0, qn, &x[0], &y[0], &z[0]))
              if(!octree->searchVector(x[nod], y[nod], z[nod], &w[0], step2,
                                       0, qn, &x[0], &y[0], &z[0]))
                octree->searchTensor(x[nod], y[nod], z[nod], &w[0], step2,
                                     0, qn, &x[0], &y[0], &z[0]);
          }
          else
            for(int comp = 0; comp < numComp2; comp++)
              data2->getValue(step2, ent, ele, nod, comp, w[comp]);
          for(int comp = 0; comp < numComp; comp++){
            if(component >= 0 && component != comp) continue;
            values[0] = x[nod]; values[1] = y[nod]; values[2] = z[nod];
            values[3] = time; values[4] = step;
            values[5] = v[comp];
            for(int i = 0; i < 9; i++) values[6 + i] = v[i];
            values[15] = w[comp];
            for(int i = 0; i < 9; i++) values[16 + i] = w[i];
            if(f.eval(values, res))
	      for(int comp = 0; comp < numComp3; comp++)
                data1->setValue(step, ent, ele, nod, comp, res[comp]);
            if(data1->isNodeData()) data1->tagNode(step, ent, ele, nod, 1);
          }
        }
      }
    }
  }

  if(octree) delete octree;

  data1->finalize();
  v1->setChanged(true);

  return v1;
}