FileInterpolatorSAFunc.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS 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.
 
   CEPS 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 CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "pde/common/functions/FileInterpolatorSAFunc.hpp"
#include "pde/common/AbstractDiscretization.hpp"
#include "pde/problem/AbstractPdeProblem.hpp"
 
#include <vtkCellLocator.h>
#include <vtkPointLocator.h>
#include <vtkUnstructuredGrid.h>
 
FileInterpolatorSAFunc::FileInterpolatorSAFunc(
  AbstractPdeProblem* pb,
  CepsReal            snapTimeRef,
  const CepsString&   filesBase,
  CepsUInt            tOrder
):
  m_problem           (pb),
  m_snapTimeRef       (snapTimeRef),
  m_tOrder            (tOrder),
  m_filesBase         (filesBase),
  m_nDigitsFile       (1),
  m_loadedTMin        (0.),
  m_loadedTMax        (0.),
  m_loadedSnapIdMin   (0),
  m_loadedSnapIdMax   (0),
  m_loadedData        ({}),
  m_loadedData0D      ({}),
  m_itpWeights        ({}),
  m_ptLocatorTolerance(1.e-8)
{
 
  CEPS_ABORT_IF(ceps::isNullPtr(pb),
    "Must pass a valid pointer to problem to create interpolator from data files"
  );
  
  // Check presence of files, and determine the number of digits in files (max is 6...)
  CepsBool foundFile = false;
  CepsUInt n         = 1;
  
  while (not foundFile and n<7)
  {
    CepsString fileName = ceps::getDir(m_filesBase) + "/snapshots/" + 
                          ceps::getBaseName(m_filesBase) + "-" +  CepsString(n,'0');
    CepsString dir      = ceps::getDir(m_filesBase);
    foundFile           = ceps::fileExists(fileName+".vtk");
    n++;
  }
 
  m_nDigitsFile = n-1;
  CEPS_ABORT_IF(not foundFile,
    "Could not find the data files with base name " << m_filesBase << 
    "***.vtk\n   from which interpolation should be made.\n" <<
    "Interpolation from files is available with vtk only."
  );
 
  m_ptLocatorTolerance = m_problem->getParameters()->getReal("convergence closest point tolerance", m_ptLocatorTolerance);
 
  computeInterpolationWeights();
 
}
 
CepsReal
FileInterpolatorSAFunc::eval(CepsStandardArgs args) 
{
  
  loadFilesForTime(args.t);
 
  CepsVector<CepsReal> itpxvals(m_loadedData.size(),0.);
  for (CepsUInt i=0; i<m_loadedData.size(); i++)
  {
    if (m_itpWeights.contains(args.dofId)) // Spatial dofs only
    {
      for (auto [vtkId,w] : m_itpWeights[args.dofId])
        itpxvals[i] += w*m_loadedData[i][args.unknownId][vtkId];
      // There might be a case, for unknowns that are not defined everywhere,
      // for which a point is located outside the subdomain of definition of the unknown, but still
      // in the mesh. In this case, the results is Nan.
      // Then we interpolate with the non nans points 
      if (std::isnan(itpxvals[i]) or std::isinf(itpxvals[i]))
      {
        CepsMap<vtkIdType,CepsReal> corrected = {};
        CepsReal                    sumW      = 0.e0;
        for (auto [vtkId,w] : m_itpWeights[args.dofId])
        {
          if (not (std::isnan(m_loadedData[i][args.unknownId][vtkId]) or 
                  (std::isinf(m_loadedData[i][args.unknownId][vtkId]))))
          {
            corrected[vtkId]  = w;
            sumW             += w;
          }
        }
        // If all nans, then nan should be interpolated, I guess... ¯\_(ツ)_/¯
        // We correct and reinterpolate only if there is an ok point
        if (corrected.size()!=0)
        {
          for (auto& [vtkId,w]: corrected)
            w /= sumW;
          m_itpWeights[args.dofId] = corrected;
          for (auto [vtkId,w] : m_itpWeights[args.dofId])
            itpxvals[i] += w*m_loadedData[i][args.unknownId][vtkId];
        }
      }
    }
    else // 0D unknowns
    {
      itpxvals[i] = m_loadedData0D[i][args.unknownId];
        ceps::checkNanOrInf(itpxvals[i]," in eval OD coeff "+ceps::toString(args.unknownId));
    }
  }
 
  return ceps::math::interpolateLagrange(args.t,itpxvals,m_loadedTMin,m_loadedTMax);
 
}
 
CepsEnum
FileInterpolatorSAFunc::getFlags() const
{
  return CepsFunctionFlag::DofIndex | CepsFunctionFlag::Time;
}
 
 
CepsEnum
ceps::getFlagsOf(FileInterpolatorSAFunc* f)
{
  return f->getFlags();
}
 
void
FileInterpolatorSAFunc::computeInterpolationWeights()
{ 
  
  auto reader = vtkSmartPointer<vtkUnstructuredGridReader>::New();  
  reader->SetFileName(getVtkFileName(0).c_str());
  reader->Update();
 
  auto pointLocator = vtkSmartPointer<vtkPointLocator>::New(); 
  auto cellLocator  = vtkSmartPointer<vtkCellLocator >::New();   
  pointLocator->SetDataSet(reader->GetOutput()); 
  cellLocator ->SetDataSet(reader->GetOutput());   
  pointLocator->BuildLocator();   
  cellLocator ->BuildLocator();  
 
  CepsReal* xRef = ceps::newArray<CepsReal>(3);
 
  AbstractDiscretization* sd = m_problem->getSpatialDiscretization();
  auto dofs = sd->getDistributedDofs();
 
  for (auto toParse: {dofs->getOwned(),dofs->getHalo()})
  for (DegreeOfFreedom* dof:toParse)
  if (ceps::isValidPtr(dof->getVertex())) // No 0D unknowns
  {
    
    CepsReal* x       = dof->getVertex()->getCoordinates().data();
    vtkIdType refPtId = pointLocator->FindClosestPoint(x);
    reader->GetOutput()->GetPoint(refPtId,xRef);
 
    CepsMap<vtkIdType,CepsReal> tmpW;
 
    // Found point coincide with parsed dof position (ideal)
    CepsReal3D xx({x[0],x[1],x[2]});
    CepsReal3D xxRef({xRef[0],xRef[1],xRef[2]});
 
    // if (ceps::approxEquals(x,xRef,3,1.e-14))
    if (areSameArray(xx,xxRef))
    {
      tmpW[refPtId] = 1.;
    }
    else // use the cell locator
    {
      double          cP[3],dist;
      int             subID,error;
      vtkIdType       refCellID;
      vtkGenericCell* g = vtkGenericCell::New();
      auto            refCellPtsIDs = vtkSmartPointer<vtkIdList>::New();
      // Get the closest cell
      cellLocator->FindClosestPoint(x,cP,g,refCellID,subID,dist);
        
      CepsVector<CepsReal> weights(g->GetNumberOfPoints());
 
      // Locate point in cell
      error = g->EvaluatePosition(x,cP,subID,xRef,dist,weights.data());
      CEPS_ABORT_IF(error ==-1 or (error==0 and dist>m_ptLocatorTolerance),
        "Error while locating points of numerical solution mesh in reference mesh cells.\n"
        "   You can increase search distance with the \"convergence closest point tolerance\"\n"
        "   parameter in the input file"
      );
 
      reader->GetOutput()->GetCellPoints(refCellID,refCellPtsIDs);        
 
      // #ifdef CEPS_DEBUG_ENABLED
      // for (vtkIdType j=0;j<g->GetNumberOfPoints();j++)
      //   ceps::checkNanOrInf(weights[j],
      //     ceps::toString(j) + "-th weight of mesh cell " + ceps::toString(refCellID) +
      //     ", when localizing dof #" + ceps::toString(dof->getGlobalIndex())
      //   );
      // #endif
 
      for (vtkIdType j=0;j<g->GetNumberOfPoints();j++)
        tmpW[refCellPtsIDs->GetId(j)] = weights[j];
      g->Delete();    
    }
 
    // Register the weights for this dof
    m_itpWeights[dof->getGlobalIndex()] = tmpW;
 
  }
 
  ceps::destroyTabular(xRef);
 
  return;
}
 
CepsString
FileInterpolatorSAFunc::getVtkFileName(CepsUInt idSnap) const
{
  std::stringstream ss;
  ss << ceps::getDir(m_filesBase) << "/snapshots/" << ceps::getBaseName(m_filesBase) 
     << "-" << std::setfill('0') << std::setw(m_nDigitsFile) << idSnap << ".vtk";
  return CepsString(ss.str());
}
 
void
FileInterpolatorSAFunc::loadFilesForTime(CepsReal t)
{
  // Determine IDs min and max for t
  CepsUInt nRefSnaps = CepsUInt(std::floor((m_problem->getEndTime()-m_problem->getStartTime())/m_snapTimeRef));
  CepsUInt idChunk   = CepsUInt(std::floor((t-m_problem->getStartTime())/(m_snapTimeRef*m_tOrder)));
 
  CepsUInt idMin     = std::max(CepsUInt(0U),std::min(idChunk*m_tOrder,nRefSnaps));
  CepsUInt idMax     = std::max(CepsUInt(0U),std::min(idMin+m_tOrder+1,nRefSnaps+1));
 
  if (idMin!=m_loadedSnapIdMin or idMax!=m_loadedSnapIdMax)
  {
    m_loadedData  .clear();
    m_loadedData0D.clear();
 
    for (CepsUInt i=idMin; i<idMax; i++)
    {
      CepsMap<CepsUnknownIndex,CepsVector<CepsReal>> tmpMap;
      CepsMap<CepsUnknownIndex,CepsReal>             tmpMap0D;
 
      auto reader = vtkSmartPointer<vtkUnstructuredGridReader>::New();  
      reader->SetFileName(getVtkFileName(i).c_str());
      reader->SetReadAllScalars(1);
      reader->SetReadAllFields(1);
      reader->Update();
 
      vtkDataArray* array = nullptr;
      for (CepsInt i=0; i<reader->GetOutput()->GetPointData()->GetNumberOfArrays(); i++)
      {
        array = reader->GetOutput()->GetPointData()->GetArray(i);
        // This check is probably useless
        if (array->GetNumberOfComponents()==1)
        {
          CepsString name = array->GetName();
          Unknown*   u    = m_problem->getUnknown(name);
          // Read only for unknowns of the problem
          if (ceps::isValidPtr(u))
          {
            vtkIdType nbData = array->GetNumberOfTuples();
            CepsVector<CepsReal> data(nbData);
            for (vtkIdType ptId = 0; ptId < nbData; ++ptId)
              array->GetTuple(ptId,&data[ptId]);
            tmpMap[u->getIdentifier()] = data;
          }
        } 
      }
      m_loadedData.push_back(tmpMap);
 
      // OD unknowns
      for (CepsInt i=0; i<reader->GetOutput()->GetFieldData()->GetNumberOfArrays(); i++)
      {
        array = reader->GetOutput()->GetFieldData()->GetArray(i);
        if (array->GetNumberOfComponents()==1 and array->GetNumberOfTuples()==1)
        {
          CepsString name = array->GetName();
          Unknown*   u    = m_problem->getUnknown(name);
          // Read only for unknowns of the problem
          if (ceps::isValidPtr(u))
          {
            CepsReal data;
            array->GetTuple(0,&data);
            tmpMap0D[u->getIdentifier()] = data;
          }
        }
      }
      m_loadedData0D.push_back(tmpMap0D);
 
    }
 
    m_loadedSnapIdMin = idMin;
    m_loadedSnapIdMax = idMax;
    m_loadedTMin = idMin*m_snapTimeRef;
    m_loadedTMax = (idMax-1)*m_snapTimeRef;
 
  }
 
  return;
}