Skip to content
Snippets Groups Projects
Select Git revision
  • fcecb101c342200368734ba25499b15bffcf4b34
  • master default protected
  • dof-renumbering
  • test-dof-hash
  • gdemesy-master-patch-30528
  • eval-space-time
  • oscillating_multiharm
  • MH_movement
  • axisqu
  • write_vtu_and_ensight_formats
  • movingband
  • CP_1972_add_vtu_file_writing
  • mortar
  • fast_freq_sweep_Resolution
  • applyresolvent_again
  • marteaua-master-patch-54323
  • patch-1
  • binde-master-patch-08072
  • binde-master-patch-52461
  • BCGSL
  • resolvent
  • getdp_3_5_0
  • getdp_3_4_0
  • getdp_3_3_0
  • getdp_3_2_0
  • getdp_3_1_0
  • getdp_3_0_4
  • getdp_3_0_3
  • getdp_3_0_2
  • getdp_3_0_1
  • getdp_3_0_0
  • onelab_mobile_2.1.0
  • getdp_2_11_3 protected
  • getdp_2_11_2 protected
  • getdp_2_11_1 protected
  • getdp_2_11_0 protected
  • getdp_2_10_0 protected
  • getdp_2_9_2 protected
  • getdp_2_9_1 protected
  • getdp_2_9_0 protected
  • getdp_2_8_0 protected
41 results

LinAlg.h

Blame
  • ModifyComponent.cpp 9.61 KiB
    // 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;
    }