Skip to content
Snippets Groups Projects
Select Git revision
  • 8c81e2edbd7d0ce6ce9966372d7cbe23992029d1
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

LineNodeBasis.h

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    NearToFarField.cpp 19.84 KiB
    // Gmsh - Copyright (C) 1997-2017 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@onelab.info>.
    //
    // Contributor(s):
    //   Ruth Sabariego
    //
    
    #include <complex>
    #include "NearToFarField.h"
    #include "OS.h"
    
    StringXNumber NearToFarFieldOptions_Number[] = {
      {GMSH_FULLRC, "Wavenumber",       NULL, 1.},
      {GMSH_FULLRC, "PhiStart",         NULL, 0.},
      {GMSH_FULLRC, "PhiEnd",           NULL, 2. * M_PI},
      {GMSH_FULLRC, "NumPointsPhi",     NULL, 60},
      {GMSH_FULLRC, "ThetaStart",       NULL, 0.},
      {GMSH_FULLRC, "ThetaEnd",         NULL, M_PI},
      {GMSH_FULLRC, "NumPointsTheta",   NULL, 30},
      {GMSH_FULLRC, "EView",            NULL, 0},
      {GMSH_FULLRC, "HView",            NULL, 1},
      {GMSH_FULLRC, "Normalize",        NULL, 1},
      {GMSH_FULLRC, "dB",               NULL, 1},
      {GMSH_FULLRC, "NegativeTime",     NULL, 0.},
      {GMSH_FULLRC, "RFar",             NULL, 0},
    };
    
    StringXString NearToFarFieldOptions_String[] = {
      {GMSH_FULLRC, "MatlabOutputFile",       NULL, "farfield.m"},
    };
    
    extern "C"
    {
      GMSH_Plugin *GMSH_RegisterNearToFarFieldPlugin()
      {
        return new GMSH_NearToFarFieldPlugin();
      }
    }
    
    std::string GMSH_NearToFarFieldPlugin::getHelp() const
    {
      return "Plugin(NearToFarField) computes the far field pattern "
        "from the near electric E and magnetic H fields on a surface "
        "enclosing the radiating device (antenna).\n\n"
        "Parameters: the wavenumber, the "
        "angular discretisation (phi in [0, 2*Pi] and theta in [0, Pi]) "
        "of the far field sphere and the indices of the views containing the "
        "complex-valued E and H fields. If `Normalize' is set, the far field "
        "is normalized to 1. If `dB' is set, the far field is computed in dB. "
        "If `NegativeTime' is set, E and H are assumed to have exp(-iwt) time "
        "dependency; otherwise they are assume to have exp(+iwt) time "
        "dependency. If `MatlabOutputFile' is given the raw far field data is "
        "also exported in Matlab format.\n\n"
        "Plugin(NearToFarField) creates one new view.";
    }
    
    int GMSH_NearToFarFieldPlugin::getNbOptions() const
    {
      return sizeof(NearToFarFieldOptions_Number) / sizeof(StringXNumber);
    }
    
    int GMSH_NearToFarFieldPlugin::getNbOptionsStr() const
    {
      return sizeof(NearToFarFieldOptions_String) / sizeof(StringXString);
    }
    
    StringXNumber *GMSH_NearToFarFieldPlugin::getOption(int iopt)
    {
      return &NearToFarFieldOptions_Number[iopt];
    }
    
    StringXString *GMSH_NearToFarFieldPlugin::getOptionStr(int iopt)
    {
      return &NearToFarFieldOptions_String[iopt];
    }
    
    // Compute field using e^{j\omega t} time dependency, following Jin in "Finite
    // Element Analysis of Antennas and Arrays", p. 176. This is not the usual `far
    // field', as it still contains the e^{ikr}/r factor.
    double GMSH_NearToFarFieldPlugin::getFarFieldJin(std::vector<element*> &allElems,
                                                     std::vector<std::vector<double> > &js,
                                                     std::vector<std::vector<double> > &ms,
                                                     double k0, double rFar, double theta,
                                                     double phi)
    {
      // theta in [0, pi] (elevation/polar angle)
      // phi in [0, 2*pi] (azimuthal angle)
    
      double sTheta = sin(theta) ; double cTheta = cos(theta) ;
      double sPhi = sin(phi) ; double cPhi = cos(phi) ;
      double r[3] = {sTheta*cPhi, sTheta*sPhi, cTheta}; // Unit vector position
    
      double Z0 = 120 * M_PI ; // free-space impedance
    
      int numComps = 3, numSteps = 2 ;
      std::vector<std::vector<double> > N  ;
      std::vector<std::vector<double> > Ns ;
      std::vector<std::vector<double> > L  ;
      std::vector<std::vector<double> > Ls ;
    
      N.resize(numSteps); Ns.resize(numSteps);
      L.resize(numSteps); Ls.resize(numSteps);
      for (int step = 0; step < numSteps; step++){
        N[step].resize(numComps, 0.); Ns[step].resize(numComps);
        L[step].resize(numComps, 0.); Ls[step].resize(numComps);
      }
    
      int i = 0 ;
      for(unsigned int ele = 0; ele < allElems.size(); ele++){
        element* e = allElems[ele] ;
        int numNodes = e->getNumNodes() ;
    
        std::vector<double > valN0(numNodes*numComps), valN1(numNodes*numComps);
        std::vector<double > valL0(numNodes*numComps), valL1(numNodes*numComps) ;
    
        for(int nod = 0; nod < numNodes; nod++){
          double x, y, z;
          e->getXYZ(nod, x, y, z) ;
          double rr, r_nod[3] = {x, y, z};
          prosca(r_nod, r, &rr) ;
          double e_jk0rr[2] = {cos(k0*rr), sin(k0*rr)} ;
    
          for(int comp = 0; comp < numComps; comp++){
            if(i < (int)js[0].size()){
              valN0[numComps * nod + comp] = js[0][i] * e_jk0rr[0] - js[1][i] * e_jk0rr[1];
              valN1[numComps * nod + comp] = js[0][i] * e_jk0rr[1] + js[1][i] * e_jk0rr[0];
              valL0[numComps * nod + comp] = ms[0][i] * e_jk0rr[0] - ms[1][i] * e_jk0rr[1];
              valL1[numComps * nod + comp] = ms[0][i] * e_jk0rr[1] + ms[1][i] * e_jk0rr[0];
              i++;
            }
          }
        }
    
        N[0][0] += e->integrate(&valN0[0], 3) ; N[1][0] += e->integrate(&valN1[0], 3) ;
        N[0][1] += e->integrate(&valN0[1], 3) ; N[1][1] += e->integrate(&valN1[1], 3) ;
        N[0][2] += e->integrate(&valN0[2], 3) ; N[1][2] += e->integrate(&valN1[2], 3) ;
    
        L[0][0] += e->integrate(&valL0[0], 3) ; L[1][0] += e->integrate(&valL1[0], 3) ;
        L[0][1] += e->integrate(&valL0[1], 3) ; L[1][1] += e->integrate(&valL1[1], 3) ;
        L[0][2] += e->integrate(&valL0[2], 3) ; L[1][2] += e->integrate(&valL1[2], 3) ;
      }
    
      // From Cartesian to spherical coordinates
      for(int step = 0; step < 2; step++){
        Ns[step][0] = N[step][0] * sTheta * cPhi + N[step][1] * sTheta * sPhi + N[step][2] *cTheta ;
        Ns[step][1] = N[step][0] * cTheta * cPhi + N[step][1] * cTheta * sPhi - N[step][2]* sTheta ;
        Ns[step][2] =-N[step][0] * sPhi          + N[step][1] * cPhi ;
    
        Ls[step][0] = L[step][0] * sTheta * cPhi + L[step][1] * sTheta * sPhi + L[step][2] *cTheta ;
        Ls[step][1] = L[step][0] * cTheta * cPhi + L[step][1] * cTheta * sPhi - L[step][2]* sTheta ;
        Ls[step][2] =-L[step][0] * sPhi          + L[step][1] * cPhi ;
      }
    
      // E_r radial component is negligible in far field
      double E_theta[2] ;
      double E_phi[2] ;
      double k0_over_4pir = k0/(4*M_PI*rFar) ;
      double cos_k0r = cos(k0*rFar) ;
      double sin_k0r = sin(k0*rFar) ;
    
      // Elevation component
      E_theta[0] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * sin_k0r -
                                     (Ls[1][2] + Z0 * Ns[1][1]) * cos_k0r ) ;
      E_theta[1] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * cos_k0r +
                                     (Ls[1][2] + Z0 * Ns[1][1]) * sin_k0r ) ;
      // Azimuthal component
      E_phi[0]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * sin_k0r -
                                     (Ls[1][1] - Z0 * Ns[1][2]) * cos_k0r ) ;
      E_phi[1]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * cos_k0r +
                                     (Ls[1][1] - Z0 * Ns[1][2]) * sin_k0r ) ;
    
      double farF =  1./2./Z0 * ( (E_theta[0] * E_theta[0] + E_theta[1] * E_theta[1]) +
                                  (E_phi[0] * E_phi[0] + E_phi[1] * E_phi[1]) ) ;
    
      return farF ;
    }
    
    // Compute far field using e^{-i\omega t} time dependency, following Monk in
    // "Finite Element Methods for Maxwell's equations", p. 233
    double GMSH_NearToFarFieldPlugin::getFarFieldMonk(std::vector<element*> &allElems,
                                                      std::vector<std::vector<double> > &ffvec,
                                                      std::vector<std::vector<double> > &js,
                                                      std::vector<std::vector<double> > &ms,
                                                      double k0, double theta, double phi)
    {
      double sTheta = sin(theta) ; double cTheta = cos(theta) ;
      double sPhi = sin(phi) ; double cPhi = cos(phi) ;
      double xHat[3] = {sTheta * cPhi, sTheta * sPhi, cTheta};
      std::complex<double> I(0., 1.);
      double Z0 = 120 * M_PI ; // free-space impedance
    
      double integral_r[3] = {0., 0., 0.}, integral_i[3] = {0., 0., 0.};
      int i = 0 ;
      for(unsigned int ele = 0; ele < allElems.size(); ele++){
        element* e = allElems[ele] ;
        int numNodes = e->getNumNodes() ;
        std::vector<double> integrand_r(numNodes * 3), integrand_i(numNodes * 3);
        for(int nod = 0; nod < numNodes; nod++){
          double y[3], xHat_dot_y;
          e->getXYZ(nod, y[0], y[1], y[2]) ;
          prosca(xHat, y, &xHat_dot_y) ;
          double n_x_e_r[3] = {-ms[0][i], -ms[0][i + 1], -ms[0][i + 2]};
          double n_x_e_i[3] = {-ms[1][i], -ms[1][i + 1], -ms[1][i + 2]};
          double n_x_h_r[3] = {js[0][i], js[0][i + 1], js[0][i + 2]};
          double n_x_h_i[3] = {js[1][i], js[1][i + 1], js[1][i + 2]};
          double n_x_h_x_xHat_r[3], n_x_h_x_xHat_i[3];
          prodve(n_x_h_r, xHat, n_x_h_x_xHat_r);
          prodve(n_x_h_i, xHat, n_x_h_x_xHat_i);
          for(int comp = 0; comp < 3; comp++){
            std::complex<double> n_x_e(n_x_e_r[comp], n_x_e_i[comp]);
            std::complex<double> n_x_h_x_xHat(n_x_h_x_xHat_r[comp], n_x_h_x_xHat_i[comp]);
            // Warning: Z0 == 1 in Monk
            std::complex<double> integrand = (n_x_e + Z0 * n_x_h_x_xHat) *
              (cos(-k0 * xHat_dot_y) + I * sin(-k0 * xHat_dot_y));
            integrand_r[3 * nod + comp] = integrand.real();
            integrand_i[3 * nod + comp] = integrand.imag();
          }
          i += 3;
        }
        for(int comp = 0; comp < 3; comp++){
          integral_r[comp] += e->integrate(&integrand_r[comp], 3);
          integral_i[comp] += e->integrate(&integrand_i[comp], 3);
        }
      }
    
      double xHat_x_integral_r[3], xHat_x_integral_i[3];
      prodve(xHat, integral_r, xHat_x_integral_r);
      prodve(xHat, integral_i, xHat_x_integral_i);
      std::complex<double> coef = I * k0 / 4. / M_PI;
      std::complex<double> einf[3] = {coef * (xHat_x_integral_r[0] + I * xHat_x_integral_i[0]),
                                      coef * (xHat_x_integral_r[1] + I * xHat_x_integral_i[1]),
                                      coef * (xHat_x_integral_r[2] + I * xHat_x_integral_i[2])};
    
      double  coef1 =  k0 / 4. / M_PI;
      for(int comp = 0; comp < 3; comp++){
        ffvec[comp][0] = -coef1 * xHat_x_integral_i[comp];
        ffvec[comp][1] =  coef1 * xHat_x_integral_r[comp];
      }
    
      return (norm(einf[0]) + norm(einf[1]) + norm(einf[2]));
    }
    
    static void printVector(FILE *fp, const std::string &name,
                            std::vector<std::vector<double> > &vec)
    {
      fprintf(fp, "%s = [", name.c_str());
      for (unsigned int i = 0; i < vec.size(); i++)
        for (unsigned int j = 0; j < vec[i].size(); j++)
          fprintf(fp, "%.16g ", vec[i][j]);
      fprintf(fp, "];\n");
    }
    
    PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
    {
      double _k0 = (double)NearToFarFieldOptions_Number[0].def;
      double _phiStart = (double)NearToFarFieldOptions_Number[1].def;
      double _phiEnd = (double)NearToFarFieldOptions_Number[2].def;
      int _NbPhi = (int)NearToFarFieldOptions_Number[3].def;
      double _thetaStart = (double)NearToFarFieldOptions_Number[4].def;
      double _thetaEnd = (double)NearToFarFieldOptions_Number[5].def;
      int _NbThe = (int)NearToFarFieldOptions_Number[6].def;
      int _eView = (int)NearToFarFieldOptions_Number[7].def;
      int _hView = (int)NearToFarFieldOptions_Number[8].def;
      bool _normalize = (bool)NearToFarFieldOptions_Number[9].def;
      bool _dB = (bool)NearToFarFieldOptions_Number[10].def;
      int _negativeTime = (int)NearToFarFieldOptions_Number[11].def;
      double _rfar = (int)NearToFarFieldOptions_Number[12].def;
    
      std::string _outFile = NearToFarFieldOptions_String[0].def;
    
      PView *ve = getView(_eView, v);
      if(!ve){
         Msg::Error("NearToFarField plugin could not find EView %i", _eView);
        return v;
      }
      PView *vh = getView(_hView, v);
      if(!vh){
        Msg::Error("NearToFarField plugin could not find HView %i", _hView);
        return v;
      }
      PViewData *eData = ve->getData() ;
      PViewData *hData = vh->getData() ;
    
      if(eData->getNumEntities() != hData->getNumEntities() ||
         eData->getNumElements() != hData->getNumElements() ||
         eData->getNumTimeSteps()!= hData->getNumTimeSteps()){
        Msg::Error("Incompatible views for E field and H field");
        return v;
      }
    
      if(eData->getNumTimeSteps() != 2 || hData->getNumTimeSteps() != 2){
        Msg::Error("Invalid number of steps for E or H fields (must be complex)");
        return v;
      }
    
      // center and radius of the visualization sphere
      SBoundingBox3d bbox = eData->getBoundingBox();
      double x0 = bbox.center().x();
      double y0 = bbox.center().y();
      double z0 = bbox.center().z();
      double lc = norm(SVector3(bbox.max(), bbox.min()));
      double r_sph = lc ? lc / 2. : 1;
    
      if(x0 != hData->getBoundingBox().center().x() ||
         y0 != hData->getBoundingBox().center().y() ||
         z0 != hData->getBoundingBox().center().z()){
        Msg::Error("E and H fields must be given on the same grid");
        return v;
      }
    
      // compute surface currents on all input elements
      std::vector<element*> allElems  ;
      std::vector<std::vector<double> > js(2) ;
      std::vector<std::vector<double> > ms(2) ;
    
      for(int ent = 0; ent < eData->getNumEntities(0); ent++){
        for(int ele = 0; ele < eData->getNumElements(0, ent); ele++){
          if(eData->skipElement(0, ent, ele)) continue;
          if(hData->skipElement(0, ent, ele)) continue;
          int numComp = eData->getNumComponents(0, ent, ele);
          if(numComp != 3) continue ;
          int dim = eData->getDimension(0, ent, ele);
          if(dim != 1 && dim != 2) continue;
          int numNodes = eData->getNumNodes(0, ent, ele);
          std::vector<double> x(numNodes), y(numNodes), z(numNodes);
          for(int nod = 0; nod < numNodes; nod++)
            eData->getNode(0, ent, ele, nod, x[nod], y[nod], z[nod]);
    
          elementFactory factory;
          allElems.push_back(factory.create(numNodes, dim, &x[0], &y[0], &z[0], true));
    
          double n[3] = {0., 0., 0.};
          if(numNodes > 2)
            normal3points(x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], n);
          else
            normal2points(x[0], y[0], z[0], x[1], y[1], z[1], n);
    
          for(int step = 0; step < 2; step++){
            for(int nod = 0; nod < numNodes; nod++){
              double h[3], e[3];
              for(int comp = 0; comp < numComp; comp++){
                eData->getValue(step, ent, ele, nod, comp, e[comp]);
                hData->getValue(step, ent, ele, nod, comp, h[comp]);
              }
              double j[3], m[3] ;
              prodve(n, h, j) ; // Js =   n x H ; Surface electric current
              prodve(e, n, m) ; // Ms = - n x E ; Surface magnetic current
              js[step].push_back(j[0]) ; js[step].push_back(j[1]) ; js[step].push_back(j[2]) ;
              ms[step].push_back(m[0]) ; ms[step].push_back(m[1]) ; ms[step].push_back(m[2]) ;
            }
          }
    
        }
      }
    
      if(allElems.empty()){
        Msg::Error("No valid elements found to compute far field");
        return v;
      }
    
      // view for far field that will contain the radiation pattern
      PView *vf = new PView();
      PViewDataList *dataFar = getDataList(vf);
    
      std::vector<std::vector<double> > phi(_NbPhi + 1), theta(_NbPhi + 1);
      std::vector<std::vector<double> > x(_NbPhi + 1), y(_NbPhi + 1), z(_NbPhi + 1);
      std::vector<std::vector<double> > farField(_NbPhi + 1);
      std::vector<std::vector<double> > farField1r(_NbPhi + 1);
      std::vector<std::vector<double> > farField2r(_NbPhi + 1);
      std::vector<std::vector<double> > farField3r(_NbPhi + 1);
      std::vector<std::vector<double> > farField1i(_NbPhi + 1);
      std::vector<std::vector<double> > farField2i(_NbPhi + 1);
      std::vector<std::vector<double> > farField3i(_NbPhi + 1);
      std::vector<std::vector<double> > farFieldVec(3);
      for(int comp = 0; comp < 3; comp++){
        farFieldVec[comp].resize(2, 0.);
      }
    
      for (int i = 0; i <= _NbPhi; i++){
        phi[i].resize(_NbThe + 1);
        theta[i].resize(_NbThe + 1);
        x[i].resize(_NbThe + 1);
        y[i].resize(_NbThe + 1);
        z[i].resize(_NbThe + 1);
        farField[i].resize(_NbThe + 1);
    
        farField1r[i].resize(_NbThe + 1);
        farField2r[i].resize(_NbThe + 1);
        farField3r[i].resize(_NbThe + 1);
        farField1i[i].resize(_NbThe + 1);
        farField2i[i].resize(_NbThe + 1);
        farField3i[i].resize(_NbThe + 1);
      }
    
      double dPhi = (_phiEnd - _phiStart) / _NbPhi;
      double dTheta = (_thetaEnd - _thetaStart) / _NbThe;
      double ffmin = 1e200, ffmax = -1e200;
      Msg::ResetProgressMeter();
      for (int i = 0; i <= _NbPhi; i++){
        for (int j = 0; j <= _NbThe; j++){
          phi[i][j] = _phiStart + i * dPhi ;
          theta[i][j] = _thetaStart + j * dTheta ;
          if(_negativeTime){
            farField[i][j] = getFarFieldMonk(allElems, farFieldVec, js, ms, _k0,
                                             theta[i][j], phi[i][j]);
            farField1r[i][j] = farFieldVec[0][0];
            farField2r[i][j] = farFieldVec[1][0];
            farField3r[i][j] = farFieldVec[2][0];
            farField1i[i][j] = farFieldVec[0][1];
            farField2i[i][j] = farFieldVec[1][1];
            farField3i[i][j] = farFieldVec[2][1];
          }
          else{
            double rfar = (_rfar ? _rfar : 10 * lc);
            farField[i][j] = getFarFieldJin(allElems, js, ms, _k0, rfar,
                                            theta[i][j], phi[i][j]);
          }
          ffmin = std::min(ffmin, farField[i][j]);
          ffmax = std::max(ffmax, farField[i][j]);
        }
        Msg::ProgressMeter(i, _NbPhi, true, "Computing far field");
      }
      for(unsigned int i = 0; i < allElems.size(); i++)
        delete allElems[i];
    
      if(_normalize){
        if(!ffmax)
          Msg::Warning("Cannot normalize far field (max = 0)");
        else
          for (int i = 0; i <= _NbPhi; i++)
            for (int j = 0; j <= _NbThe; j++)
              farField[i][j] /= ffmax ;
      }
    
      if(_dB){
        ffmin = 1e200;
        ffmax = -1e200;
        for (int i = 0; i <= _NbPhi; i++){
          for (int j = 0; j <= _NbThe; j++){
            farField[i][j] = 10 * log10(farField[i][j]);
            ffmin = std::min(ffmin, farField[i][j]);
            ffmax = std::max(ffmax, farField[i][j]);
          }
        }
      }
    
      for (int i = 0; i <= _NbPhi; i++){
        for (int j = 0; j <= _NbThe; j++){
          double df = (ffmax - ffmin);
          if(!df){
            Msg::Warning("zero far field range");
            df = 1.;
          }
          double f = (farField[i][j] - ffmin) / df; // in [0,1]
          x[i][j] = x0 + r_sph * f * sin(theta[i][j]) * cos(phi[i][j]);
          y[i][j] = y0 + r_sph * f * sin(theta[i][j]) * sin(phi[i][j]);
          z[i][j] = z0 + r_sph * f * cos(theta[i][j]);
        }
      }
    
      if(_outFile.size()){
        FILE *fp = Fopen(_outFile.c_str(), "w");
        if(fp){
          printVector(fp, "phi", phi);
          printVector(fp, "theta", theta);
          printVector(fp, "farField", farField);
    
          if(_negativeTime){
            printVector(fp, "farField1r", farField1r);
            printVector(fp, "farField2r", farField2r);
            printVector(fp, "farField3r", farField3r);
            printVector(fp, "farField1i", farField1i);
            printVector(fp, "farField2i", farField2i);
            printVector(fp, "farField3i", farField3i);
          }
    
          printVector(fp, "x", x);
          printVector(fp, "y", y);
          printVector(fp, "z", z);
          fclose(fp);
        }
        else
          Msg::Error("Could not open file '%s'", _outFile.c_str());
      }
    
      for (int i = 0; i < _NbPhi; i++){
        for (int j = 0; j < _NbThe; j++){
          if(_NbPhi == 1 || _NbThe == 1){
            dataFar->NbSP++;
            dataFar->SP.push_back(x[i][j]); dataFar->SP.push_back(y[i][j]);
            dataFar->SP.push_back(z[i][j]); dataFar->SP.push_back(farField[i][j]);
          }
          else{
            double P1[3] = {x[i][j], y[i][j], z[i][j]};
            double P2[3] = {x[i+1][j], y[i+1][j], z[i+1][j]};
            double P3[3] = {x[i+1][j+1], y[i+1][j+1], z[i+1][j+1]};
            double P4[3] = {x[i][j+1], y[i][j+1], z[i][j+1]};
            dataFar->NbSQ++;
            dataFar->SQ.push_back(P1[0]); dataFar->SQ.push_back(P2[0]);
            dataFar->SQ.push_back(P3[0]); dataFar->SQ.push_back(P4[0]);
            dataFar->SQ.push_back(P1[1]); dataFar->SQ.push_back(P2[1]);
            dataFar->SQ.push_back(P3[1]); dataFar->SQ.push_back(P4[1]);
            dataFar->SQ.push_back(P1[2]); dataFar->SQ.push_back(P2[2]);
            dataFar->SQ.push_back(P3[2]); dataFar->SQ.push_back(P4[2]);
            dataFar->SQ.push_back(farField[i][j]);
            dataFar->SQ.push_back(farField[i+1][j]);
            dataFar->SQ.push_back(farField[i+1][j+1]);
            dataFar->SQ.push_back(farField[i][j+1]);
          }
        }
      }
    
      dataFar->setName("_NearToFarField");
      dataFar->setFileName("_NearToFarField.pos");
      dataFar->finalize();
    
      return vf;
    }