Skip to content
Snippets Groups Projects
Select Git revision
  • cf473f274db07e19eec25f8c008243cf25173de8
  • master default protected
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

MeshSubEntities.h

Blame
  • F_Gmsh.cpp 7.32 KiB
    // GetDP - Copyright (C) 1997-2018 P. Dular and C. Geuzaine, University of Liege
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <getdp@onelab.info>.
    
    #include "GetDPConfig.h"
    #include "ProData.h"
    #include "F.h"
    #include "Message.h"
    
    extern struct CurrentData Current ;
    
    #if defined(HAVE_GMSH)
    
    #include <gmsh/GmshGlobal.h>
    #include <gmsh/PView.h>
    #include <gmsh/PViewData.h>
    
    void F_Field(F_ARG)
    {
      if(A->Type != VECTOR){
        Message::Error("Field[] expects XYZ coordinates as argument");
        return;
      }
      if(PView::list.empty()){
        Message::Error("No views available to interpolate from");
        return;
      }
      double x = A->Val[0];
      double y = A->Val[1];
      double z = A->Val[2];
    
      if(Fct->NbrArguments > 1){
        Message::Error("Time and additional arguments are not supported in Field: "
                       "use {Scalar,Vector,Tensor}Field instead");
        return;
      }
    
      for (int k = 0; k < Current.NbrHar; k++)
        for (int j = 0; j < 9; j++)
          V->Val[MAX_DIM * k + j] = 0. ;
      V->Type = SCALAR;
    
      std::vector<int> iview;
      if(!Fct->NbrParameters){
        // use last view by default
        iview.push_back(PView::list.back()->getTag());
      }
      else{
        for(int i = 0; i < Fct->NbrParameters; i++)
          iview.push_back(Fct->Para[i]);
      }
    
      double N = 0.;
    
      // add the values from all specified views
      for(unsigned int i = 0; i < iview.size(); i++){
    
        PView *v = PView::getViewByTag(iview[i]);
        if(!v){
          Message::Error("View with tag %d does not exist", iview[i]);
          return;
        }
    
        PViewData *data = v->getData();
    
        std::vector<double> val(9 * data->getNumTimeSteps());
        if(data->searchScalar(x, y, z, &val[0])){
          V->Val[0] += val[0];
          if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1)
            V->Val[MAX_DIM] += val[1];
          V->Type = SCALAR;
          N += 1.;
        }
        else if(data->searchVector(x, y, z, &val[0])){
          for(int j = 0; j < 3; j++)
            V->Val[j] += val[j];
          if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){
            for(int j = 0; j < 3; j++)
              V->Val[MAX_DIM + j] += val[3 + j];
          }
          V->Type = VECTOR;
          N += 1.;
        }
        else if(data->searchTensor(x, y, z, &val[0])){
          for(int j = 0; j < 9; j++)
            V->Val[j] += val[j];
          if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){
            for(int j = 0; j < 9; j++)
              V->Val[MAX_DIM + j] += val[9 + j];
          }
          V->Type = TENSOR;
          N += 1.;
        }
        else{
          Message::Error("Did not find data at point (%g,%g,%g) in View with tag %d",
                         x, y, z, iview[i]);
        }
      }
    
      if(N > 1.){
        Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
        for (int k = 0; k < Current.NbrHar; k++)
          for (int j = 0; j < 9; j++)
            V->Val[MAX_DIM * k + j] = 0. ;
      }
    }
    
    static void F_X_Field(F_ARG, int type, bool complex, bool grad=false)
    {
      if(A->Type != VECTOR){
        Message::Error("Field[] expects XYZ coordinates as argument");
        return;
      }
      if(PView::list.empty()){
        Message::Error("No views available to interpolate from");
        return;
      }
    
      double x = A->Val[0];
      double y = A->Val[1];
      double z = A->Val[2];
      int numComp = (type == SCALAR) ? (grad ? 3 : 1) :
        (type == VECTOR) ? (grad ? 9 : 3) : 9; // TODO: grad of tensor
    
      int NbrArg = Fct->NbrArguments ;
      int TimeStep = 0, MatchElement = 0;
      if(NbrArg >= 2){
        if((A+1)->Type != SCALAR){
          Message::Error("Expected scalar second argument (time step)");
          return;
        }
        TimeStep = (int)(A+1)->Val[0];
      }
      if(NbrArg >= 3){
        if((A+2)->Type != SCALAR){
          Message::Error("Expected scalar second argument (element matching flag)");
          return;
        }
        MatchElement = (int)(A+2)->Val[0];
      }
    
      // TODO: we could treat the third arguement as a tolerance (and call
      // searchScalarWithTol & friends)
    
      // Complex{Scalar,Vector,Tensor}Field assume that the Gmsh view contains real
      // and imaginary parts for each step
      if(complex) TimeStep *= 2;
    
      for (int k = 0; k < Current.NbrHar; k++)
        for (int j = 0; j < numComp; j++)
          V->Val[MAX_DIM * k + j] = 0. ;
      V->Type = (numComp == 1) ? SCALAR : (numComp == 3) ? VECTOR : TENSOR;
    
      std::vector<int> iview;
      if(!Fct->NbrParameters){
        // use last view by default
        iview.push_back(PView::list.back()->getTag());
      }
      else{
        for(int i = 0; i < Fct->NbrParameters; i++)
          iview.push_back(Fct->Para[i]);
      }
    
      int qn = 0;
      double *qx = 0, *qy = 0, *qz = 0;
      if(Current.Element){
        qn = MatchElement ? Current.Element->GeoElement->NbrNodes : 0;
        qx = Current.Element->x;
        qy = Current.Element->y;
        qz = Current.Element->z;
      }
    
      double N = 0.;
    
      // add the values from all specified views
      for(unsigned int i = 0; i < iview.size(); i++){
        PView *v = PView::getViewByTag(iview[i]);
        if(!v){
          Message::Error("View with tag %d does not exist", iview[i]);
          return;
        }
        PViewData *data = v->getData();
        if(TimeStep < 0 || TimeStep >= data->getNumTimeSteps()){
          Message::Error("Invalid step %d in View with tag %d", TimeStep, iview[i]);
          continue;
        }
        std::vector<double> val(numComp * data->getNumTimeSteps());
        bool found = false;
        switch(type){
        case SCALAR :
          if(data->searchScalar(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
            found = true;
          break;
        case VECTOR :
          if(data->searchVector(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
            found = true;
          break;
        case TENSOR :
          // TODO: grad of tensor not allowed yet - not sure how we should return
          // the values; provide 3 routines that return 3 tensors, or add argumemt
          // to select what to return?
          if(data->searchTensor(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, false))
            found = true;
          break;
        }
        if(found){
          for(int j = 0; j < numComp; j++)
            V->Val[j] += val[numComp * TimeStep + j];
          if(complex && Current.NbrHar == 2 && data->getNumTimeSteps() > TimeStep + 1)
            for(int j = 0; j < numComp; j++)
              V->Val[MAX_DIM + j] += val[numComp * (TimeStep + 1) + j];
          N += 1.;
        }
      }
    
      if(N > 1.){
        Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
        for (int k = 0; k < Current.NbrHar; k++)
          for (int j = 0; j < numComp; j++)
            V->Val[MAX_DIM * k + j] /= N ;
      }
    }
    
    #else
    
    void F_Field(F_ARG)
    {
      Message::Error("You need to compile GetDP with Gmsh support to use 'Field'");
      V->Val[0] = 0. ;
      V->Type = SCALAR ;
    }
    
    static void F_X_Field(F_ARG, int type, bool complex, bool grad=false)
    {
      Message::Error("You need to compile GetDP with Gmsh support to use 'Field'");
      V->Val[0] = 0. ;
      V->Type = SCALAR ;
    }
    
    #endif
    
    void F_ScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false); }
    void F_VectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false); }
    void F_TensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, false); }
    void F_ComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true); }
    void F_ComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true); }
    void F_ComplexTensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, true); }
    
    void F_GradScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false, true); }
    void F_GradVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false, true); }
    void F_GradComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true, true); }
    void F_GradComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true, true); }