MeasureVP2D.cpp

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: MeasureVP2D.cpp                                                  |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D is free software: you can redistribute it and/or modify it             |
| under the terms of the GNU General Public License as published by           |
| the Free Software Foundation, either version 3 of the License, or           |
| (at your option) any later version.                                         |
|                                                                             |
| VM2D is distributed in the hope that it will be useful, but WITHOUT         |
| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       |
| FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License       |
| for more details.                                                           |
|                                                                             |
| You should have received a copy of the GNU General Public License           |
| along with VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
 
#if defined(_WIN32)
#include <direct.h>
#endif
 
#include "MeasureVP2D.h"
 
#include "Airfoil2D.h"
#include "Boundary2D.h"
 
#include "Mechanics2D.h"
#include "Mechanics2DDeformable.h"
#include "Preprocessor.h"
#include "StreamParser.h"
#include "Velocity2D.h"
#include "Wake2D.h"
#include "World2D.h"
#include "Gmres2D.h"
 
using namespace VM2D;
 
//Конструктор
MeasureVP::MeasureVP(const World2D& W_)
    : W(W_)
{
    sstrCounter = 0;
    wakeVP.reset(new WakeDataBase(W_));
};//MeasureVP(...)
 
//Чтение точек, в которых нужно посчитать давление и скорость
void MeasureVP::ReadPointsFromFile(const std::string& dir)
{
    std::string filename = dir + "pointsVP";
    std::ifstream VPFile;
 
    if (fileExistTest(filename, W.getInfo(), true, { "txt", "TXT" }))
    {
        std::stringstream VPFile(VMlib::Preprocessor(filename).resultString);
 
        VMlib::StreamParser VPParser(W.getInfo(), "velocity & pressure parser", VPFile);
        
        VPParser.get("points", initialPoints);
        VPParser.get("history", historyPoints);
 
        for (auto& pt : initialPoints)
            pt = pt.rotated(-W.getPassport().rotateAngleVpPoints * PI / 180.0);
        
        for (auto& pt : historyPoints)
            pt = pt.rotated(-W.getPassport().rotateAngleVpPoints * PI / 180.0);
 
 
        if (historyPoints.size() > 0)
        {
            sstr.resize(historyPoints.size());
            sstrCounter = 0;
 
            VMlib::CreateUserDirectory(dir, "velPres");
            std::string VPFileNameList = W.getPassport().dir + "velPres/listPoints";
            std::ofstream VPFileList(VPFileNameList.c_str());
 
            VMlib::PrintLogoToTextFile(VPFileList, VPFileNameList.c_str(), "List of points, where velocity and pressure are measured in csv-files");
 
            VMlib::PrintHeaderToTextFile(VPFileList, "No.     FileName      X     Y");
 
            for (size_t q = 0; q < historyPoints.size(); ++q)
            {
                VPFileList << '\n' << q << '\t' << VMlib::fileNameStep("VP-atPoint-", 2, q, "csv") << '\t' << historyPoints[q][0] << '\t' << historyPoints[q][1];
 
                std::string VPFileNameCsv;
                VPFileNameCsv = W.getPassport().dir + "velPres/" + VMlib::fileNameStep("VP-atPoint-", 2, q, "csv");
 
                std::ofstream VPFileCsv(VPFileNameCsv.c_str());
 
                if (!W.getPassport().calcCoefficients)
                    VPFileCsv << "pt,t,Vx,Vy,p" << std::endl;
                else
                    VPFileCsv << "pt,t,CVx,CVy,Cp" << std::endl;
                
                VPFileCsv.close();
                VPFileCsv.clear();
            }
 
            VPFileList.close();
            VPFileList.clear();
        }
    }
}//ReadPointsFromFile(...)
 
//Чтение точек, в которых нужно посчитать давление и скорость
void MeasureVP::SetPoints(const std::vector<Point2D>& points, const std::vector<Point2D>& history)
{
    initialPoints = points;
    historyPoints = history;
    
    sstr.resize(historyPoints.size());
    sstrCounter = 0;
 
    for (auto& pt : initialPoints)
        pt = pt.rotated(-W.getPassport().rotateAngleVpPoints * PI / 180.0);
 
    for (auto& pt : historyPoints)
        pt = pt.rotated(-W.getPassport().rotateAngleVpPoints * PI / 180.0);
 
 
    if (historyPoints.size() > 0)
    {
        VMlib::CreateUserDirectory(W.getPassport().dir, "velPres");
        std::string VPFileNameList = W.getPassport().dir + "velPres/listPoints";
        std::ofstream VPFileList(VPFileNameList.c_str());
 
        VMlib::PrintLogoToTextFile(VPFileList, VPFileNameList.c_str(), "List of points, where velocity and pressure are measured in csv-files");
 
        VMlib::PrintHeaderToTextFile(VPFileList, "No.     FileName      X     Y");
 
        for (size_t q = 0; q < historyPoints.size(); ++q)
        {
            VPFileList << '\n' << q << '\t' << VMlib::fileNameStep("VP-atPoint-", 2, q, "csv") << '\t' << historyPoints[q][0] << '\t' << historyPoints[q][1];
 
            std::string VPFileNameCsv;
            VPFileNameCsv = W.getPassport().dir + "velPres/" + VMlib::fileNameStep("VP-atPoint-", 2, q, "csv");
 
            std::ofstream VPFileCsv(VPFileNameCsv.c_str());
 
            if (!W.getPassport().calcCoefficients)
                VPFileCsv << "pt,t,Vx,Vy,p" << std::endl;
            else
                VPFileCsv << "pt,t,CVx,CVy,Cp" << std::endl;
 
            VPFileCsv.close();
            VPFileCsv.clear();
        }
 
        VPFileList.close();
        VPFileList.clear();
    }
}//SetPoints(...)
 
 
// Инициализация векторов для вычисления скоростей и давлений
void MeasureVP::Initialization()
{
    W.getTimers().start("VelPres");
 
    if ((W.getPassport().timeDiscretizationProperties.saveVPstep > 0) && (!(W.getCurrentStep() % W.getPassport().timeDiscretizationProperties.saveVPstep)))
    {
        //W.getInfo('i') << "Preparing VP points" << std::endl;
 
        wakeVP->vtx.clear();
        wakeVP->vtx.reserve(initialPoints.size() + historyPoints.size() + elasticPoints.size());
 
        Vortex2D addvtx;
 
        velocity.clear();
        pressure.clear();
        domainRadius.clear();
 
        //отключаем проверку того, что точки расположены вне профиля
        /*
        //определяем точки, которые не находятся внутри профилей
        bool inside;
        for (size_t i = 0; i < initialPoints.size(); ++i)
        {
            inside = false;
            for (size_t j = 0; j < W.getNumberOfAirfoil(); ++j)
            {
                if (W.getAirfoil(j).isInsideGabarits(initialPoints[i]) && W.getAirfoil(j).IsPointInAirfoil(initialPoints[i]))
                {
                    inside = true;
                    break;
                }
            }
            if (!inside)
            {
                addvtx.r() = initialPoints[i];
                wakeVP.vtx.push_back(addvtx);
            }
        }
        */
 
        //добавляем все точки, а не только те, которые вне профиля
        for (size_t i = 0; i < initialPoints.size(); ++i)
        {
            addvtx.r() = initialPoints[i];
            addvtx.g() = 0.0;
            addvtx.sigma() = -1.0;
            wakeVP->vtx.push_back(addvtx);
            
        }//for i
 
        for (size_t i = 0; i < historyPoints.size(); ++i)
        {
            addvtx.r() = historyPoints[i];
            addvtx.g() = 0.0;
            addvtx.sigma() = -1.0;
            wakeVP->vtx.push_back(addvtx);
        }//for i
 
 
        elasticPoints.clear();
                
        for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
        {
            MechanicsDeformable* ptr = dynamic_cast<MechanicsDeformable*>(&(W.getNonConstMechanics(afl)));
            if ((ptr != nullptr) && (ptr->beam->fsi))
            {
                if (afl > 0)
                {
                    W.getInfo('e') << "Only airfoil 0 can be elastic!" << std::endl;
                    exit(-345);
                }               
 
                //if (ptr->beam->fsi) //turek
                {
                    elasticPoints.reserve(ptr->chord.size() * 2);
                    for (size_t q = 0; q < ptr->chord.size(); ++q)
                    {
                        std::pair<size_t, size_t> pnls = ptr->chord[q].infPanels;
                        elasticPoints.push_back(
                            0.5 * (W.getAirfoil(afl).getR(pnls.first) + W.getAirfoil(afl).getR(pnls.first + 1))
                            + 0.01 * W.getAirfoil(afl).len[pnls.first] * W.getAirfoil(afl).nrm[pnls.first]
                        );
 
                        elasticPoints.push_back(
                            0.5 * (W.getAirfoil(afl).getR(pnls.second) + W.getAirfoil(afl).getR(pnls.second + 1))
                            + 0.01 * W.getAirfoil(afl).len[pnls.second] * W.getAirfoil(afl).nrm[pnls.second]
                        );
                    }
 
                    //std::ofstream elastPointsPositionsFile(W.getPassport().dir + "elastPointsPos-" + std::to_string(W.getCurrentStep()) + ".txt");
                    //for (size_t q = 0; q < elasticPoints.size(); ++q)
                    //  elastPointsPositionsFile << elasticPoints[q][0] << " " << elasticPoints[q][1] << std::endl;
                    //elastPointsPositionsFile.close();
                }                           
            }
            else //fish and other types
            {
                elasticPoints.reserve(W.getAirfoil(afl).getNumberOfPanels());
                for (size_t q = 0; q < W.getAirfoil(afl).getNumberOfPanels(); ++q)
                    elasticPoints.push_back(0.5 * (W.getAirfoil(afl).getR(q) + W.getAirfoil(afl).getR(q + 1)) + W.getAirfoil(afl).nrm[q] * 1e-10);
            }
        }
 
        /*
        std::cout << "Elast_points added!\n";
        for (size_t afl = 0; afl < W.getNumberOfAirfoil(); ++afl)
        {
            elasticPoints.reserve(W.getAirfoil(afl).getNumberOfPanels());
            for (size_t q = 0; q < W.getAirfoil(afl).getNumberOfPanels(); ++q)
                elasticPoints.push_back(0.5 * (W.getAirfoil(afl).getR(q) + W.getAirfoil(afl).getR(q + 1)) + W.getAirfoil(afl).nrm[q] * 1e-4);
        }
        */
 
        for (size_t i = 0; i < elasticPoints.size(); ++i)
        {
            addvtx.r() = elasticPoints[i];
            addvtx.g() = 0.0;
            addvtx.sigma() = -1.0;
            wakeVP->vtx.push_back(addvtx);
        }//for i
 
    }
    W.getTimers().stop("VelPres");
 
}//Initialization()
 
//Расчет поля давления
void MeasureVP::CalcPressure()
{
    W.getTimers().start("VelPres");
 
    int addWSize = (int)wakeVP->vtx.size(); 
    pressure.resize(addWSize);  
    
    Point2D dri;
    Point2D Vi;
    Point2D vi;
 
    const Point2D& V0 = W.getV0();
    const double& dt = W.getPassport().timeDiscretizationProperties.dt;
 
    double P0 = /*1*/ 0.5 * V0.length2();
 
    Point2D cPan;
 
#pragma warning (push)
#pragma warning (disable: 4101)
    double alpha;
    double lambda;
    double dst2eps;
#pragma warning (pop)
 
#pragma omp parallel for default(none) private(alpha, lambda, dri, Vi, vi, dst2eps, cPan) shared(P0, dt, V0, addWSize, std::cout, IDPI) 
    for (int i = 0; i < addWSize; ++i)
    {
        const Point2D& pt = wakeVP->vtx[i].r();
                
        pressure[i] = P0;
        pressure[i] -= 0.5 * velocity[i].length2(); //2
 
        
        for (size_t j = 0; j < W.getWake().vtx.size(); ++j)
        {
            dri = pt - W.getWake().vtx[j].r();
            Vi = W.getVelocity().wakeVortexesParams.convVelo[j] + W.getVelocity().wakeVortexesParams.diffVelo[j] + V0;
 
            const double sigma2 = sqr(W.getWake().vtx[j].sigma());
            dst2eps = VMlib::boundDenom(dri.length2(), sigma2);
 
            vi = IDPI * W.getWake().vtx[j].g() / dst2eps * dri.kcross();
            pressure[i] += vi & Vi; //3
        }
 
 
        for (size_t q = 0; q < W.getNumberOfAirfoil(); ++q)
        //size_t q = 0;
        {
            const Airfoil& afl = W.getAirfoil(q);
            const Boundary& bou = W.getBoundary(q);                     
 
            for (size_t pnl = 0; pnl < afl.getNumberOfPanels(); ++pnl)
            {
                Point2D s = pt - afl.getR(pnl);
                Point2D p = pt - afl.getR(pnl + 1);
                Point2D u0 = afl.tau[pnl];
                double alpha = VMlib::Alpha(s, p);
                double lambda = VMlib::Lambda(s, p);
 
                Point2D skos = IDPI * ((alpha * u0).kcross() + (lambda * u0));
 
                Vi = 0.5 * (bou.sheets.freeVortexSheet(pnl, 0) /* - afl.gammaThrough[pnl] / afl.len[pnl]*/) * u0;
                vi = skos.kcross() * bou.sheets.freeVortexSheet(pnl, 0);
                //if ((i == 0) && (pnl == 0))
                //  std::cout << "vi & Vi = " << (vi & Vi) << std::endl;
 
                pressure[i] += vi & Vi; //4
 
                Vi = 0.5 * (afl.getV(pnl) + afl.getV(pnl + 1));
                vi = skos.kcross() * (bou.sheets.attachedVortexSheet(pnl, 0) - bou.oldSheets.attachedVortexSheet(pnl, 0));
                pressure[i] += vi & Vi; //4a
 
                vi = skos * (bou.sheets.attachedSourceSheet(pnl, 0) - bou.oldSheets.attachedSourceSheet(pnl, 0));
                pressure[i] += vi & Vi; //4b
            }
        }
 
        for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
        {
            const Point2D& rcm = W.getAirfoil(bou).rcm;
 
 
            //
            //if ((i==0) && (bou == 2))
            //{
            //  for (int q = 0; q < W.getAirfoil(bou).possibleWays.size(); ++q)
            //  {
            //      std::cout << "Possible way #" << q << ": " << std::endl;
            //      for (size_t s = 0; s < W.getAirfoil(bou).possibleWays[q].size(); ++s)
            //          std::cout << W.getAirfoil(bou).possibleWays[q][s] << " ";
            //      std::cout << std::endl;
            //  }
            //  
            //  std::ofstream of(W.getPassport().dir + "ways2.txt");
            //  for (size_t q = 0; q < W.getAirfoil(bou).getNumberOfPanels(); ++q)
            //  {
            //      cPan = 0.5 * (W.getAirfoil(bou).getR(q) + W.getAirfoil(bou).getR(q + 1));
 
            //      alpha = 0.0;
 
            //      int way = W.getAirfoil(bou).wayToVertex[q];
            //      if (way == 0)
            //          of << q << " " << rcm[0] << " " << rcm[1] << " " << cPan[0] << " " << cPan[1] << std::endl;
            //      else
            //      {
            //          of << q << " " << rcm[0] << " " << rcm[1];
            //          for (int q = 0; q < W.getAirfoil(bou).possibleWays[way - 1].size() + 1; ++q)
            //          {
            //              Point2D start = ((q == 0) ? rcm : W.getAirfoil(bou).possibleWays[way - 1][q - 1]);
            //              Point2D finish = ((q == W.getAirfoil(bou).possibleWays[way - 1].size()) ? cPan : W.getAirfoil(bou).possibleWays[way - 1][q]);
            //              of << " " << finish[0] << " " << finish[1];
            //          }
            //          of << std::endl;
            //      }
            //  }               
            //  of.close();
            //}
            //
 
 
            for (size_t j = 0; j < W.getAirfoil(bou).getNumberOfPanels(); ++j)
            {
                cPan = 0.5 * (W.getAirfoil(bou).getR(j) + W.getAirfoil(bou).getR(j + 1));
                
                alpha = 0.0;
                
                int way = W.getAirfoil(bou).wayToVertex[j];
                //if (way < 0)              
                //  W.getInfo('i') << "Pressure computation is incorrect, way inside airfoil is undefined" << std::endl;                
 
                if (way == 0)
                    alpha = atan2((cPan - pt) ^ (rcm - pt), (cPan - pt) & (rcm - pt));
                else
                    if (way > 0)
                    {
                        for (int q = 0; q < W.getAirfoil(bou).possibleWays[way - 1].size() + 1; ++q)
                        {
                            Point2D start = ((q == 0) ? rcm : W.getAirfoil(bou).possibleWays[way - 1][q - 1]);
                            Point2D finish = ((q == W.getAirfoil(bou).possibleWays[way - 1].size()) ? cPan : W.getAirfoil(bou).possibleWays[way - 1][q]);
                            alpha += atan2((finish - pt) ^ (start - pt), (finish - pt) & (start - pt));
                        }
                    }
                
 
                pressure[i] += IDPI * alpha * (W.getBoundary(bou).sheets.freeVortexSheet(j, 0) + W.getBoundary(bou).sheets.attachedVortexSheet(j, 0) - W.getBoundary(bou).oldSheets.attachedVortexSheet(j, 0)) * W.getAirfoil(bou).len[j] / dt;
 
                pressure[i] -= IDPI * alpha * W.getAirfoil(bou).gammaThrough[j] / dt;
 
                lambda = 0.5 * log((pt - cPan).length2());
                pressure[i] -= IDPI * lambda * (W.getBoundary(bou).sheets.attachedSourceSheet(j, 0) - W.getBoundary(bou).oldSheets.attachedSourceSheet(j, 0)) * W.getAirfoil(bou).len[j] / dt;
            }
        }//for bou      
    }//for i
    
    for (size_t i = 0; i < pressure.size(); ++i)
        pressure[i] *= W.getPassport().physicalProperties.rho;
 
    W.getTimers().stop("VelPres");
}//CalcPressure()
 
 
//Расчет поля давления
#ifdef USE_CUDA
void MeasureVP::GPUCalcPressure()
{
    W.getTimers().start("VelPres");
 
    int addWSize = (int)wakeVP->vtx.size();
    const Point2D& V0 = W.getV0();
 
#if (defined(USE_CUDA))
    {
        cuCopyFixedArray(W.getWake().devVelPtr, W.getNonConstVelocity().wakeVortexesParams.convVelo.data(), 2 * sizeof(double) * W.getWake().vtx.size());
        //W.getNonConstCuda().RefreshWake();
    }
#endif
 
    size_t npnl = 0;
    for (size_t q = 0; q < W.getNumberOfAirfoil(); ++q)
        npnl += W.getAirfoil(q).getNumberOfPanels();
 
    std::vector<double> oldAttVortexSheet, oldAttSourceSheet, gammaThroughIntensity;
    oldAttVortexSheet.reserve(npnl);
    oldAttSourceSheet.reserve(npnl);
    gammaThroughIntensity.reserve(npnl);
    
    for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
        for (size_t p = 0; p < W.getBoundary(bou).afl.getNumberOfPanels(); ++p)
        {
            oldAttVortexSheet.push_back(W.getBoundary(bou).oldSheets.attachedVortexSheet(p, 0));
            oldAttSourceSheet.push_back(W.getBoundary(bou).oldSheets.attachedSourceSheet(p, 0));
            gammaThroughIntensity.push_back(W.getAirfoil(bou).gammaThrough[p] / W.getAirfoil(bou).len[p]);
        }
 
    double* devOldVortexSheet;
    cudaMalloc(&devOldVortexSheet, npnl * sizeof(double));
    cudaMemcpy(devOldVortexSheet, oldAttVortexSheet.data(), npnl * sizeof(double), cudaMemcpyHostToDevice);
 
    double* devOldSourceSheet;
    cudaMalloc(&devOldSourceSheet, npnl * sizeof(double));
    cudaMemcpy(devOldSourceSheet, oldAttSourceSheet.data(), npnl * sizeof(double), cudaMemcpyHostToDevice);
 
    double* devGammaThroughIntensity;
    cudaMalloc(&devGammaThroughIntensity, npnl * sizeof(double));
    cudaMemcpy(devGammaThroughIntensity, gammaThroughIntensity.data(), npnl * sizeof(double), cudaMemcpyHostToDevice);
    
    double* devDiffVelo;
    cudaMalloc(&devDiffVelo, W.getWake().vtx.size() * 2 * sizeof(double));
    cudaMemcpy(devDiffVelo, W.getVelocity().wakeVortexesParams.diffVelo.data(), W.getWake().vtx.size() * 2 * sizeof(double), cudaMemcpyHostToDevice);
 
 
    cuCalculatePressure(addWSize, wakeVP->devVtxPtr, 
        W.getWake().vtx.size(), W.getWake().devVtxPtr, V0[0], V0[1], W.getWake().devVelPtr, devDiffVelo,
        npnl, W.getAirfoil(0).devRPtr, W.getAirfoil(0).devFreeVortexSheetPtr, W.getAirfoil(0).devAttachedVortexSheetPtr, W.getAirfoil(0).devAttachedSourceSheetPtr,
        devOldVortexSheet, devOldSourceSheet, devGammaThroughIntensity,
        this->devPressurePtr);
    
    cudaFree(devOldVortexSheet);
    cudaFree(devOldSourceSheet);
    cudaFree(devGammaThroughIntensity);
    cudaFree(devDiffVelo);
 
    std::vector<double> presGPU(addWSize);
    cuCopyMemFromDev(presGPU.data(), this->devPressurePtr, addWSize * sizeof(double), 18);
 
    //printf("-------------------\n");
 
    pressure.resize(addWSize);
 
    Point2D dri;
    Point2D Vi;
    Point2D vi;
 
    
    const double& dt = W.getPassport().timeDiscretizationProperties.dt;
 
    double P0 = /*1*/ 0.5 * V0.length2();
 
    Point2D cPan;
 
#pragma warning (push)
#pragma warning (disable: 4101)
    double alpha;
    double lambda;
    double dst2eps;
#pragma warning (pop)
 
    for (int i = 0; i < addWSize; ++i)
    {
        //if (fabs(pressure[i] - presGPU[i]) > 1e-5)
        //  printf("different pressure: %d, %f, %f\n", i, pressure[i], presGPU[i]);
        pressure[i] = presGPU[i];
    }
 
#pragma omp parallel for default(none) private(alpha, lambda, dri, Vi, vi, dst2eps, cPan) shared(P0, dt, V0, addWSize, std::cout, IDPI) 
    for (int i = 0; i < addWSize; ++i)
    {
        const Point2D& pt = wakeVP->vtx[i].r();
        for (size_t bou = 0; bou < W.getNumberOfBoundary(); ++bou)
        {
            const Point2D& rcm = W.getAirfoil(bou).rcm;
 
 
            /*
            if ((i==0) && (bou == 2))
            {
                for (int q = 0; q < W.getAirfoil(bou).possibleWays.size(); ++q)
                {
                    std::cout << "Possible way #" << q << ": " << std::endl;
                    for (size_t s = 0; s < W.getAirfoil(bou).possibleWays[q].size(); ++s)
                        std::cout << W.getAirfoil(bou).possibleWays[q][s] << " ";
                    std::cout << std::endl;
                }
 
                std::ofstream of(W.getPassport().dir + "ways2.txt");
                for (size_t q = 0; q < W.getAirfoil(bou).getNumberOfPanels(); ++q)
                {
                    cPan = 0.5 * (W.getAirfoil(bou).getR(q) + W.getAirfoil(bou).getR(q + 1));
 
                    alpha = 0.0;
 
                    int way = W.getAirfoil(bou).wayToVertex[q];
                    if (way == 0)
                        of << q << " " << rcm[0] << " " << rcm[1] << " " << cPan[0] << " " << cPan[1] << std::endl;
                    else
                    {
                        of << q << " " << rcm[0] << " " << rcm[1];
                        for (int q = 0; q < W.getAirfoil(bou).possibleWays[way - 1].size() + 1; ++q)
                        {
                            Point2D start = ((q == 0) ? rcm : W.getAirfoil(bou).possibleWays[way - 1][q - 1]);
                            Point2D finish = ((q == W.getAirfoil(bou).possibleWays[way - 1].size()) ? cPan : W.getAirfoil(bou).possibleWays[way - 1][q]);
                            of << " " << finish[0] << " " << finish[1];
                        }
                        of << std::endl;
                    }
                }
                of.close();
            }
            */
 
 
            for (size_t j = 0; j < W.getAirfoil(bou).getNumberOfPanels(); ++j)
            {
                cPan = 0.5 * (W.getAirfoil(bou).getR(j) + W.getAirfoil(bou).getR(j + 1));
 
                alpha = 0.0;
 
                int way = W.getAirfoil(bou).wayToVertex[j];
                //if (way < 0)              
                //  W.getInfo('i') << "Pressure computation is incorrect, way inside airfoil is undefined" << std::endl;                
 
                if (way == 0)
                    alpha = atan2((cPan - pt) ^ (rcm - pt), (cPan - pt) & (rcm - pt));
                else
                    if (way > 0)
                    {
                        for (int q = 0; q < W.getAirfoil(bou).possibleWays[way - 1].size() + 1; ++q)
                        {
                            Point2D start = ((q == 0) ? rcm : W.getAirfoil(bou).possibleWays[way - 1][q - 1]);
                            Point2D finish = ((q == W.getAirfoil(bou).possibleWays[way - 1].size()) ? cPan : W.getAirfoil(bou).possibleWays[way - 1][q]);
                            alpha += atan2((finish - pt) ^ (start - pt), (finish - pt) & (start - pt));
                        }
                    }
 
                pressure[i] += IDPI * alpha * (W.getBoundary(bou).sheets.freeVortexSheet(j, 0) + W.getBoundary(bou).sheets.attachedVortexSheet(j, 0) - W.getBoundary(bou).oldSheets.attachedVortexSheet(j, 0)) * W.getAirfoil(bou).len[j] / dt;
 
                pressure[i] -= IDPI * alpha * W.getAirfoil(bou).gammaThrough[j] / dt;
 
                lambda = 0.5 * log((pt - cPan).length2());
                pressure[i] -= IDPI * lambda * (W.getBoundary(bou).sheets.attachedSourceSheet(j, 0) - W.getBoundary(bou).oldSheets.attachedSourceSheet(j, 0)) * W.getAirfoil(bou).len[j] / dt;
            }
        }//for bou
    }//for i
 
    for (size_t i = 0; i < pressure.size(); ++i)
    {
        pressure[i] += P0 - 0.5 * velocity[i].length2(); //2;
        pressure[i] *= W.getPassport().physicalProperties.rho;
    }
 
    W.getTimers().stop("VelPres");
}//GPUCalcPressure()
#endif
 
//Сохранение в файл вычисленных скоростей и давлений
void MeasureVP::SaveVP()
{
    W.getTimers().start("Save");
 
    double scaleV = 1.0, scaleP = 1.0;
    if (W.getPassport().calcCoefficients)
    {
        const double& vRef = W.getPassport().physicalProperties.vRef;
        scaleV = 1.0 / vRef;
        scaleP = 1.0 / (0.5 * W.getPassport().physicalProperties.rho * sqr(vRef));
    }
 
 
    std::ofstream outfile;
 
    VMlib::CreateUserDirectory(W.getPassport().dir, "velPres");
    const size_t nRealPressurePoints = initialPoints.size() + historyPoints.size(); //число реальных точек - без учета тех, что для гидроупругости
 
    if (W.getPassport().timeDiscretizationProperties.fileTypeVP.second == 0) //text format VTK
    {
        std::string fname = VMlib::fileNameStep("VelPres", W.getPassport().timeDiscretizationProperties.nameLength, W.getCurrentStep(), "vtk");
        outfile.open(W.getPassport().dir + "velPres/" + fname);
 
        outfile << "# vtk DataFile Version 2.0\n";
        outfile << "VM2D VTK result: " << (W.getPassport().dir + "velPres/" + fname).c_str() << " saved " << VMlib::CurrentDataTime() << '\n';
        outfile << "ASCII\n";
        outfile << "DATASET UNSTRUCTURED_GRID\n";
        outfile << "POINTS " << nRealPressurePoints << " float\n";
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            double xi = (wakeVP->vtx[i].r())[0];
            double yi = (wakeVP->vtx[i].r())[1];
            outfile << xi << " " << yi << " " << "0.0\n";
        }//for i
 
        outfile << "CELLS " << nRealPressurePoints << " " << 2 * nRealPressurePoints << '\n';
        for (size_t i = 0; i < nRealPressurePoints; ++i)
            outfile << "1 " << i << '\n';
 
        outfile << "CELL_TYPES " << nRealPressurePoints << '\n';
        for (size_t i = 0; i < nRealPressurePoints; ++i)
            outfile << "1\n";
 
        outfile << '\n';
        outfile << "POINT_DATA " << nRealPressurePoints << '\n';
 
        if (!W.getPassport().calcCoefficients)
            outfile << "VECTORS V float\n";
        else
            outfile << "VECTORS CV float\n";
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            outfile << velocity[i][0] * scaleV << " " << velocity[i][1] * scaleV << " 0.0\n";
        }//for i
 
        outfile << '\n';
 
        if (!W.getPassport().calcCoefficients)
            outfile << "SCALARS P float 1\n";
        else
            outfile << "SCALARS CP float 1\n";
 
        outfile << "LOOKUP_TABLE default\n";
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            outfile << pressure[i] * scaleP << '\n';
        }//for i
 
        outfile.close();
    }
    else if (W.getPassport().timeDiscretizationProperties.fileTypeVP.second == 1) //binary format VTK
    {
        //Тест способа хранения чисел
        uint16_t x = 0x0001;
        bool littleEndian = (*((uint8_t*)&x));
        const char eolnBIN[] = "\n";
 
        std::string fname = VMlib::fileNameStep("VelPres", W.getPassport().timeDiscretizationProperties.nameLength, W.getCurrentStep(), "vtk");
        outfile.open(W.getPassport().dir + "velPres/" + fname, std::ios::out | std::ios::binary);
 
        outfile << "# vtk DataFile Version 3.0" << "\r\n" << "VM2D VTK result: " << (W.getPassport().dir + "velPres/" + fname).c_str() << " saved " << VMlib::CurrentDataTime() << eolnBIN;
        outfile << "BINARY" << eolnBIN;
        outfile << "DATASET UNSTRUCTURED_GRID" << eolnBIN << "POINTS " << nRealPressurePoints << " " << "float" << eolnBIN;
 
 
        Eigen::VectorXf pData = Eigen::VectorXf::Zero(nRealPressurePoints);
        Eigen::VectorXf vData = Eigen::VectorXf::Zero(nRealPressurePoints * 3);
        Eigen::VectorXf rData = Eigen::VectorXf::Zero(nRealPressurePoints * 3);
 
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            rData(3 * i + 0) = (float)(wakeVP->vtx[i].r())[0];
            rData(3 * i + 1) = (float)(wakeVP->vtx[i].r())[1];
        }//for i
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            vData(3 * i + 0) = (float)(velocity[i][0] * scaleV);
            vData(3 * i + 1) = (float)(velocity[i][1] * scaleV);
        }//for i
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            pData(i) = (float)(pressure[i] * scaleP);
        }//for i
 
        //POINTS
        if (littleEndian)
            for (int i = 0; i < nRealPressurePoints * 3; ++i)
                VMlib::SwapEnd(rData(i));
        outfile.write(reinterpret_cast<char*>(rData.data()), nRealPressurePoints * 3 * sizeof(float));
 
        // CELLS
        std::vector<int> cells(2 * nRealPressurePoints);
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            cells[2 * i + 0] = 1;
            cells[2 * i + 1] = (int)i;
        }
 
        std::vector<int> cellsTypes;
        cellsTypes.resize(nRealPressurePoints, 1);
 
        if (littleEndian)
        {
            for (int i = 0; i < nRealPressurePoints * 2; ++i)
                VMlib::SwapEnd(cells[i]);
 
            for (int i = 0; i < nRealPressurePoints; ++i)
                VMlib::SwapEnd(cellsTypes[i]);
        }
 
        outfile << eolnBIN << "CELLS " << nRealPressurePoints << " " << nRealPressurePoints * 2 << eolnBIN;
        outfile.write(reinterpret_cast<char*>(cells.data()), nRealPressurePoints * 2 * sizeof(int));
        outfile << eolnBIN << "CELL_TYPES " << nRealPressurePoints << eolnBIN;
        outfile.write(reinterpret_cast<char*>(cellsTypes.data()), nRealPressurePoints * sizeof(int));
 
        //VECTORS V
        if (littleEndian)
            for (int i = 0; i < nRealPressurePoints * 3; ++i)
                VMlib::SwapEnd(vData(i));
 
        outfile << eolnBIN << "POINT_DATA " << nRealPressurePoints << eolnBIN;
 
        if (!W.getPassport().calcCoefficients)
            outfile << "VECTORS V " << "float" << eolnBIN;
        else
            outfile << "VECTORS CV " << "float" << eolnBIN;
 
        outfile.write(reinterpret_cast<char*>(vData.data()), nRealPressurePoints * 3 * sizeof(float));
 
        //SCALARS P 
        if (littleEndian)
            for (int i = 0; i < nRealPressurePoints; ++i)
                VMlib::SwapEnd(pData(i));
 
        if (!W.getPassport().calcCoefficients)
            outfile << eolnBIN << "SCALARS P " << "float" << " 1" << eolnBIN;
        else
            outfile << eolnBIN << "SCALARS CP " << "float" << " 1" << eolnBIN;
 
        outfile << "LOOKUP_TABLE default" << eolnBIN;
        outfile.write(reinterpret_cast<char*>(pData.data()), nRealPressurePoints * sizeof(float));
 
        outfile << eolnBIN;
 
        outfile.close();
    }
    else if (W.getPassport().timeDiscretizationProperties.fileTypeVP.second == 2) //csv
    {
        std::string fname = VMlib::fileNameStep("VelPres", W.getPassport().timeDiscretizationProperties.nameLength, W.getCurrentStep(), "csv");
        outfile.open(W.getPassport().dir + "velPres/" + fname);
 
        if (!W.getPassport().calcCoefficients)
            outfile << "point,x,y,Vx,Vy,P\n";
        else
            outfile << "point,x,y,CVx,CVy,CP\n";
 
        for (size_t i = 0; i < nRealPressurePoints; ++i)
        {
            double xi = (wakeVP->vtx[i].r())[0];
            double yi = (wakeVP->vtx[i].r())[1];
            outfile << i << "," << xi << "," << yi \
                << "," << velocity[i][0] * scaleV << "," << velocity[i][1] * scaleV \
                << "," << pressure[i] * scaleP << '\n';
        }//for i
        outfile.close();
    }
    else if (W.getPassport().timeDiscretizationProperties.fileTypeVP.second == 3) //csvBundle
    {
        std::string fnameBunCsv = W.getPassport().dir + "velPres/" + "velPresBundle.csv";
        std::ofstream VPFileBunCsv(fnameBunCsv.c_str(), W.getCurrentStep() ? std::ios::app : std::ios::out);
        {
            if (W.getCurrentStep() == 0)
                VPFileBunCsv << nRealPressurePoints << '\n';
 
 
            for (size_t q = 0; q < nRealPressurePoints; ++q)
            {
                VPFileBunCsv << W.getCurrentStep() << "," \
                    << W.getCurrentTime() << "," \
                    << q << "," \
                    << wakeVP->vtx[q].r()[0] << "," << wakeVP->vtx[q].r()[1] << "," \
                    << velocity[q][0] * scaleV << "," << velocity[q][1] * scaleV << "," \
                    << pressure[q] * scaleP << '\n';
            }
        }
 
    }
 
 
    //Вывод в csv-файлы точек "historyPoints"
    #pragma omp parallel for
    for (int q = 0; q < (int)historyPoints.size(); ++q)
    {
 
        if (W.getCurrentStep() == 0)
        {
            std::string VPFileNameCsv;
            VPFileNameCsv = W.getPassport().dir + "velPres/" + VMlib::fileNameStep("VP-atPoint-", 2, q, "csv");
 
            std::ofstream VPFileCsv(VPFileNameCsv.c_str(), W.getCurrentStep() ? std::ios::app : std::ios::out);
 
            if (!W.getPassport().calcCoefficients)
                VPFileCsv << "point,time,Vx,Vy,P\n";
            else
                VPFileCsv << "point,time,CVx,CVy,CP\n";
 
            VPFileCsv.close();
        }
 
        sstr[q] << q << "," << W.getCurrentTime() << ","
            << velocity[initialPoints.size() + q][0] * scaleV << ","
            << velocity[initialPoints.size() + q][1] * scaleV << ","
            << pressure[initialPoints.size() + q] * scaleP << '\n';
    
        if (q == 0)
            ++sstrCounter;
    }
        
    if (sstrCounter == 100)
    {
        #pragma omp parallel for
        for (int q = 0; q < (int)historyPoints.size(); ++q)
        {
            std::string VPFileNameCsv;
            VPFileNameCsv = W.getPassport().dir + "velPres/" + VMlib::fileNameStep("VP-atPoint-", 2, q, "csv");
 
            std::ofstream VPFileCsv(VPFileNameCsv.c_str(), W.getCurrentStep() ? std::ios::app : std::ios::out);
 
            VPFileCsv << sstr[q].str();
 
            VPFileCsv.close();
 
            sstr[q].str("");
            sstr[q].clear();
 
            
            if (q == 0)
                sstrCounter = 0;
        }
    }
        
    W.getTimers().stop("Save");
}//SaveVP()
 
std::vector<std::pair<Point2D, double>> MeasureVP::GetVPinElasticPoints()
{
    std::vector<std::pair<Point2D, double>> result;
    result.reserve(elasticPoints.size());
    
    const size_t nRealPressurePoints = initialPoints.size() + historyPoints.size(); //число реальных точек - без учета тех, что для гидроупругости
 
    for (size_t q = 0; q < elasticPoints.size(); ++q)
        result.push_back({ velocity[nRealPressurePoints + q], pressure[nRealPressurePoints + q] });
    
    return result;
}