JPartDataBi4.cpp

//HEAD_DSCODES
/*
 <DUALSPHYSICS>  Copyright (c) 2020 by Dr Jose M. Dominguez et al. (see http://dual.sphysics.org/index.php/developers/). 

 EPHYSLAB Environmental Physics Laboratory, Universidade de Vigo, Ourense, Spain.
 School of Mechanical, Aerospace and Civil Engineering, University of Manchester, Manchester, U.K.

 This file is part of DualSPHysics. 

 DualSPHysics is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License 
 as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version.
 
 DualSPHysics 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 Lesser General Public License for more details. 

 You should have received a copy of the GNU Lesser General Public License along with DualSPHysics. If not, see <http://www.gnu.org/licenses/>. 
*/

/// \file JPartDataBi4.cpp \brief Implements the class \ref JPartDataBi4.

#include "JPartDataBi4.h"
//#include "JBinaryData.h"
#include "JPartDataHead.h"
#include "Functions.h"
#include <fstream>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <sstream>

using namespace std;

//##############################################################################
//# JPartDataBi4
//##############################################################################
//==============================================================================
/// Constructor.
//==============================================================================
JPartDataBi4::JPartDataBi4(){
  ClassName="JPartDataBi4";
  Data=NULL;
  Reset();
}

//==============================================================================
/// Destructor.
//==============================================================================
JPartDataBi4::~JPartDataBi4(){
  DestructorActive=true;
  delete Data; Data=NULL;
}

//==============================================================================
/// Initialisation of variables.
//==============================================================================
void JPartDataBi4::Reset(){
  ResetData();
  Dir="";
  Piece=0;
  Npiece=1;
}

//==============================================================================
/// Elimina informacion de Data.
/// Deletes information from data.
//==============================================================================
void JPartDataBi4::ResetData(){
  delete Data; 
  Data=new JBinaryData(ClassName);
  Part=Data->CreateItem("Part");
  Cpart=0;
}

//==============================================================================
/// Elimina informacion de PARTs.
/// Deletes information from PARTs.
//==============================================================================
void JPartDataBi4::ResetPart(){
  Part->Clear();
}

//==============================================================================
/// Devuelve la memoria reservada.
/// Returns allocated memory
//==============================================================================
long long JPartDataBi4::GetAllocMemory()const{  
  long long s=0;
  s+=Data->GetAllocMemory();
  return(s);
}

//==============================================================================
/// Devuelve nombre de fichero PART segun los parametros indicados.
/// Returns the filename PART according to the specified parameters.
//==============================================================================
std::string JPartDataBi4::GetFileNamePart(unsigned cpart,unsigned piece,unsigned npiece){
  string fname="Part";
  if(npiece>1)fname=fname+fun::PrintStr("_p%02d",piece);
  return(fname+fun::PrintStr("_%04u.bi4",cpart));
}

//==============================================================================
/// Devuelve nombre de fichero de caso segun los parametros indicados.
/// Returns filename's case according to the specified parameters.
//==============================================================================
std::string JPartDataBi4::GetFileNameCase(const std::string &casename,unsigned piece,unsigned npiece){
  string fname=casename;
  if(npiece>1)fname=fname+fun::PrintStr("_p%02d",piece);
  return(fname+".bi4");
}

//==============================================================================
/// Devuelve nombre de fichero info de caso segun los parametros indicados.
/// Returns filename's case info according to the specified parameters.
//==============================================================================
std::string JPartDataBi4::GetFileNameInfo(unsigned piece,unsigned npiece){
  string fname="PartInfo";
  if(npiece>1)fname=fname+fun::PrintStr("_p%02d",piece);
  return(fname+".ibi4");
}

//==============================================================================
/// Devuelve nombre de fichero encontrado segun configuracion e indica si esta o 
/// no dividido en varias piezas (0:No se encontro, 1:Una pieza, 2:Varias piezas).
/// Returns file name found depending on configuration and indicates 
/// whether or not it is divided into several parts
/// (0:not found, 1:a piece, 2:Several parts).
//==============================================================================
std::string JPartDataBi4::GetFileData(std::string casename,std::string dirname,unsigned cpart,byte &npiece){
  byte npie=0;
  string file;
  if(casename.empty()){
    dirname=fun::GetDirWithSlash(dirname);
    if(fun::FileExists(dirname+JPartDataBi4::GetFileNamePart(cpart,0,1))){
      file=dirname+JPartDataBi4::GetFileNamePart(cpart,0,1);
      npie=1;
    }
    else if(fun::FileExists(dirname+JPartDataBi4::GetFileNamePart(cpart,0,2))){
      file=dirname+JPartDataBi4::GetFileNamePart(cpart,0,2);
      npie=2;
    }
  }
  else{
    if(fun::FileExists(JPartDataBi4::GetFileNameCase(casename,0,1))){
      file=JPartDataBi4::GetFileNameCase(casename,0,1);
      npie=1;
    }
    else if(fun::FileExists(JPartDataBi4::GetFileNameCase(casename,0,2))){
      file=JPartDataBi4::GetFileNameCase(casename,0,2);
      npie=2;
    }
  }
  npiece=npie;
  return(file);
}

//==============================================================================
/// Object configuration from JPartDataHead object.
//==============================================================================
void JPartDataBi4::Config(unsigned piece,unsigned npiece,std::string dir,const JPartDataHead* phead){
  ConfigBasic(piece,npiece,phead->GetRunCode(),phead->GetAppName(),phead->GetCaseName()
    ,phead->GetData2d(),phead->GetData2dPosY(),dir);
  ConfigParticles(phead->GetCaseNp(),phead->GetCaseNfixed(),phead->GetCaseNmoving()
    ,phead->GetCaseNfloat(),phead->GetCaseNfluid(),phead->GetCasePosMin()
    ,phead->GetCasePosMax(),phead->GetNpDynamic(),phead->GetReuseIds());
  ConfigCtes(phead->GetDp(),phead->GetH(),phead->GetB(),phead->GetRhopZero()
    ,phead->GetGamma(),phead->GetMassBound(),phead->GetMassFluid());
  ConfigSimMap(phead->GetMapPosMin(),phead->GetMapPosMax());
  ConfigSimPeri(phead->GetPeriMode(),phead->GetPeriXinc(),phead->GetPeriYinc(),phead->GetPeriZinc());
  ConfigSymmetry(phead->GetSymmetry());
  ConfigSplitting(phead->GetSplitting());
}

//==============================================================================
/// Configuracion de variables basicas.
/// Configuration of basic variables.
//==============================================================================
void JPartDataBi4::ConfigBasic(unsigned piece,unsigned npiece,std::string runcode
  ,std::string appname,std::string casename,bool data2d,double data2dposy,const std::string &dir)
{
  ResetData();
  Piece=piece; Npiece=npiece;
  Dir=fun::GetDirWithSlash(dir);
  Data->SetvUint("Piece",Piece);
  Data->SetvUint("Npiece",Npiece);
  Data->SetvText("RunCode",runcode);
  Data->SetvText("Date",fun::GetDateTime());
  Data->SetvText("AppName",appname);
  Data->SetvText("CaseName",casename);
  Data->SetvBool("Data2d",data2d);
  Data->SetvDouble("Data2dPosY",data2dposy);
  ConfigSimMap(TDouble3(0),TDouble3(0));
  ConfigSimPeri(PERI_Unknown,TDouble3(0),TDouble3(0),TDouble3(0));
  ConfigSimDiv(DIV_Unknown);
}

//==============================================================================
/// Configuracion de numero de particulas y dominio del caso.
/// Setting number of particles and domain of the case.
//==============================================================================
void JPartDataBi4::ConfigParticles(ullong casenp,ullong casenfixed,ullong casenmoving
  ,ullong casenfloat,ullong casenfluid,tdouble3 caseposmin,tdouble3 caseposmax
  ,bool npdynamic,bool reuseids)
{
  if(casenp!=casenfixed+casenmoving+casenfloat+casenfluid)Run_Exceptioon("Error in the number of particles.");
  Data->SetvUllong("CaseNp",casenp);
  Data->SetvUllong("CaseNfixed",casenfixed);
  Data->SetvUllong("CaseNmoving",casenmoving);
  Data->SetvUllong("CaseNfloat",casenfloat);
  Data->SetvUllong("CaseNfluid",casenfluid);
  Data->SetvDouble3("CasePosMin",caseposmin);
  Data->SetvDouble3("CasePosMax",caseposmax);
  Data->SetvBool("NpDynamic",npdynamic);
  Data->SetvBool("ReuseIds",reuseids);
}

//==============================================================================
/// Configuracion de constantes.
/// Configuration of constants.
//==============================================================================
void JPartDataBi4::ConfigCtes(double dp,double h,double b,double rhop0,double gamma
  ,double massbound,double massfluid)
{
  Data->SetvDouble("Dp",dp);
  Data->SetvDouble("H",h);
  Data->SetvDouble("B",b);
  Data->SetvDouble("Rhop0",rhop0);
  Data->SetvDouble("Gamma",gamma);
  Data->SetvDouble("MassBound",massbound);
  Data->SetvDouble("MassFluid",massfluid);
}

//==============================================================================
/// Configuracion de variables de simulacion: map limits.
/// Configuration of variables of simulation: map limits.
//==============================================================================
void JPartDataBi4::ConfigSimMap(tdouble3 mapposmin,tdouble3 mapposmax){
  Data->SetvDouble3("MapPosMin",mapposmin);
  Data->SetvDouble3("MapPosMax",mapposmax);
}

//==============================================================================
/// Configuracion de variables de condiciones periodicas.
/// Configuration of variables of periodic conditions.
//==============================================================================
void JPartDataBi4::ConfigSimPeri(TpPeri tperi,tdouble3 perixinc,tdouble3 periyinc,tdouble3 perizinc){
  Data->SetvInt("PeriMode",int(tperi));
  Data->SetvDouble3("PeriXinc",perixinc);
  Data->SetvDouble3("PeriYinc",periyinc);
  Data->SetvDouble3("PeriZinc",perizinc);
}

//==============================================================================
/// Devuelve configuracion de ejes periodicos.
/// Returns configuration of periodic axes.
//==============================================================================
TpPeri JPartDataBi4::Get_PeriMode()const{
  string varname="PeriMode";
  //-Maintains backwards compatibility with previous formats (27-04-2018).
  if(!GetData()->ExistsValue("PeriMode",JBinaryDataDef::DatInt) && GetData()->ExistsValue("PeriActive",JBinaryDataDef::DatInt))varname="PeriActive";
  return((TpPeri)GetData()->GetvInt(varname)); 
}

//==============================================================================
/// Configuracion de variables de simulacion: axis division.
/// Configuration of variables of simulation: axis division.
//==============================================================================
void JPartDataBi4::ConfigSimDiv(TpAxisDiv axisdiv){
  Data->SetvInt("AxisDiv",int(axisdiv));
}

//==============================================================================
/// Configuracion de variables de simetria con respecto al plano y=0.
/// Configuration of variables of symmetry according plane y=0.
//==============================================================================
void JPartDataBi4::ConfigSymmetry(bool symmetry){
  Data->SetvBool("Symmetry",symmetry);
}

//==============================================================================
/// Configuracion uso de Splitting.
/// Configuration used for Splitting.
//==============================================================================
void JPartDataBi4::ConfigSplitting(bool splitting){
  Data->SetvBool("Splitting",splitting);
}

//==============================================================================
/// Devuelve nombre de part segun su numero.
/// Returns name of part according to their number.
//==============================================================================
std::string JPartDataBi4::GetNamePart(unsigned cpart){
  return(fun::PrintStr("PART_%04u",cpart));
}

//==============================================================================
/// Anhade informacion de nuevo part.
// Add information to new part.
//==============================================================================
JBinaryData* JPartDataBi4::AddPartInfo(unsigned cpart,double timestep,unsigned npok,unsigned nout,unsigned step,double runtime,tdouble3 domainmin,tdouble3 domainmax,ullong nptotal,ullong idmax){
  Part->Clear();
  Cpart=cpart;
  Part->SetName(GetNamePart(cpart));
  Part->SetvUint("Cpart",cpart);
  Part->SetvDouble("TimeStep",timestep);
  Part->SetvUint("Npok",npok);
  Part->SetvUint("Nout",nout);
  Part->SetvUint("Step",step);
  Part->SetvDouble("RunTime",runtime);
  Part->SetvDouble3("DomainMin",domainmin);
  Part->SetvDouble3("DomainMax",domainmax);
  if(nptotal)Part->SetvUllong("NpTotal",nptotal);
  if(idmax)Part->SetvUllong("IdMax",idmax);
  return(Part);
}

//==============================================================================
/// Anhade datos (definidos por el usuario) de particulas de de nuevo part.
/// Add data (defined by user) of particles to new part.
//==============================================================================
void JPartDataBi4::AddPartDataVar(const std::string &name,JBinaryDataDef::TpData type,unsigned npok,const void *v,bool externalpointer){
  if(!v)Run_Exceptioon("The pointer data is invalid.");
  //-Comprueba valor de npok. Checks value of npok.
  if(Part->GetvUint("Npok")!=npok)Run_Exceptioon("Part information is invalid.");
  //-Crea array con particulas validas. Creates valid particles array.
  Part->CreateArray(name,type,npok,v,externalpointer);
}

//==============================================================================
/// Anhade datos (definidos por el usuario) de particulas de de nuevo part.
/// Add data (defined by user) of particles to new part.
//==============================================================================
void JPartDataBi4::AddPartData(const std::string &name,unsigned npok,const void *v
  ,TpTypeData type,bool externalpointer)
{
  switch(type){
    case TypeUchar:    AddPartData(name,npok,(byte    *)v,externalpointer);   break;
    case TypeUshort:   AddPartData(name,npok,(word    *)v,externalpointer);   break;
    case TypeUint:     AddPartData(name,npok,(unsigned*)v,externalpointer);   break;
    case TypeFloat:    AddPartData(name,npok,(float   *)v,externalpointer);   break;
    case TypeDouble:   AddPartData(name,npok,(double  *)v,externalpointer);   break;
    case TypeUint3:    AddPartData(name,npok,(tuint3  *)v,externalpointer);   break;
    case TypeFloat3:   AddPartData(name,npok,(tfloat3 *)v,externalpointer);   break;
    case TypeDouble3:  AddPartData(name,npok,(tdouble3*)v,externalpointer);   break;
    default: Run_Exceptioon("Type of pointer is unknown.");
  }
}

//==============================================================================
/// Anhade datos de particulas de de nuevo part.
/// Adds data of particles to new part.
//==============================================================================
void JPartDataBi4::AddPartData(unsigned npok,const unsigned *idp,const ullong *idpd
  ,const tfloat3 *pos,const tdouble3 *posd,const tfloat3 *vel,const float *rhop
  ,bool externalpointer)
{
  if(!idp&&!idpd)Run_Exceptioon("The id of particles is invalid.");
  if(!pos&&!posd)Run_Exceptioon("The position of particles is invalid.");
  //-Comprueba valor de npok. Checks value of npok.
  if(Part->GetvUint("Npok")!=npok)Run_Exceptioon("Part information is invalid.");
  //-Crea array con particulas validas. Creates valid particles array.
  if(idpd)Part->CreateArray("Idpd",JBinaryDataDef::DatUllong,npok,idpd,externalpointer);
  else    Part->CreateArray("Idp" ,JBinaryDataDef::DatUint,npok,idp,externalpointer);
  if(posd)Part->CreateArray("Posd",JBinaryDataDef::DatDouble3,npok,posd,externalpointer);
  else    Part->CreateArray("Pos" ,JBinaryDataDef::DatFloat3,npok,pos,externalpointer);
  Part->CreateArray("Vel",JBinaryDataDef::DatFloat3,npok,vel,externalpointer);
  Part->CreateArray("Rhop",JBinaryDataDef::DatFloat,npok,rhop,externalpointer);
}

//==============================================================================
/// Anhade datos Splitting de particulas de de nuevo part.
/// Add data Splitting of particles to new part.
//==============================================================================
void JPartDataBi4::AddPartDataSplitting(unsigned npok,const float *mass
  ,const float *hvar,bool externalpointer)
{
  if(!mass || !hvar)Run_Exceptioon("The pointer data is invalid.");
  //-Comprueba valor de npok. Checks value of npok.
  if(Part->GetvUint("Npok")!=npok)Run_Exceptioon("Part information is invalid.");
  if(!Data->GetvBool("Splitting"))Run_Exceptioon("Splitting is not configured.");
  //-Crea array con particulas validas. Creates valid particles array.
  Part->CreateArray("Mass",JBinaryDataDef::DatFloat,npok,mass,externalpointer);
  Part->CreateArray("Hvar",JBinaryDataDef::DatFloat,npok,hvar,externalpointer);
}

//==============================================================================
/// Graba le fichero BI4 indicado.
/// Writes indicated BI4 file.
//==============================================================================
void JPartDataBi4::SaveFileData(std::string fname){
  //-Comprueba que Part tenga algun array de datos. Check that Part has array with data.
  if(!Part->GetArraysCount())Run_Exceptioon("There is not array of particles data.");
  //-Graba fichero. Record file.
  Data->SaveFile(Dir+fname,false,true);
  Part->RemoveArrays();
  //Data->SaveFileXml(Dir+fun::GetWithoutExtension(fname)+"__.xml");
}

//==============================================================================
/// Graba fichero BI4 con el nombre da caso indicado.
/// Writes file BI4 with the case name indicated.
//==============================================================================
void JPartDataBi4::SaveFileCase(std::string casename){
  SaveFileData(GetFileNameCase(casename,Piece,Npiece));
}

//==============================================================================
/// Graba fichero PART con datos de particulas.
/// Writes file PART with data of particles.
//==============================================================================
void JPartDataBi4::SaveFilePart(){
  SaveFileData(GetFileNamePart(Cpart,Piece,Npiece));
}

//==============================================================================
/// Graba info de PART.
/// Writes info PART.
//==============================================================================
void JPartDataBi4::SaveFileInfo(){
  Data->SetHideItems(true,false);
  Part->SetHideArrays(true,false);
  Part->SaveFileListApp(Dir+GetFileNameInfo(Piece,Npiece),ClassName+"_Info",true,false);
  Data->SetHideItems(false,false);
}

//==============================================================================
/// Devuelve el numero de piezas del fichero indicado.
/// Returns the number of parts from the indicated file.
//==============================================================================
unsigned JPartDataBi4::GetPiecesFile(std::string file)const{
  unsigned npieces=0;
  if(fun::FileExists(file)){
    JBinaryData dat(ClassName);
    dat.OpenFileStructure(file,ClassName);
    npieces=dat.GetvUint("Npiece");
  }
  return(npieces);
}

//==============================================================================
/// Devuelve el numero de piezas del caso indicado.
/// Returns the number of parts of the case indicated.
//==============================================================================
unsigned JPartDataBi4::GetPiecesFileCase(std::string dir,std::string casename)const{
  unsigned npieces=0;
  if(fun::FileExists(dir+GetFileNameCase(casename,0,1)))npieces=1;
  else npieces=GetPiecesFile(dir+GetFileNameCase(casename,0,2));
  return(npieces);
}

//==============================================================================
/// Devuelve el numero de piezas del caso indicado.
/// Returns the number of parts of the case indicated.
//==============================================================================
unsigned JPartDataBi4::GetPiecesFilePart(std::string dir,unsigned cpart)const{
  unsigned npieces=0;
  if(fun::FileExists(Dir+GetFileNamePart(cpart,0,1)))npieces=1;
  else npieces=GetPiecesFile(dir+GetFileNamePart(cpart,0,2));
  return(npieces);
}

//==============================================================================
/// Graba fichero BI4 con el nombre da caso indicado.
/// Writes file BI4 with the case name indicated.
//==============================================================================
void JPartDataBi4::LoadFileData(std::string file,unsigned cpart,unsigned piece,unsigned npiece){
  ResetData();
  Cpart=cpart; Piece=piece; Npiece=npiece;
  Data->OpenFileStructure(file,ClassName);
  if(Piece!=Data->GetvUint("Piece")||Npiece!=Data->GetvUint("Npiece"))Run_Exceptioon("PART configuration is invalid.");
  Part=Data->GetItem(GetNamePart(Cpart));
  if(!Part)Run_Exceptioon("PART data is invalid.");
  Cpart=Part->GetvUint("Cpart");
}

//==============================================================================
/// Carga fichero BI4 con el nombre da caso indicado.
/// Load file BI4 with the case name indicated.
//==============================================================================
void JPartDataBi4::LoadFileCase(std::string dir,std::string casename,unsigned piece,unsigned npiece){
  LoadFileData(fun::GetDirWithSlash(dir)+GetFileNameCase(casename,piece,npiece),0,piece,npiece);
}

//==============================================================================
/// Carga fichero PART con datos de particulas.
/// Load file PART with data of particles.
//==============================================================================
void JPartDataBi4::LoadFilePart(std::string dir,unsigned cpart,unsigned piece,unsigned npiece){
  LoadFileData(fun::GetDirWithSlash(dir)+GetFileNamePart(cpart,piece,npiece),cpart,piece,npiece);
}

//==============================================================================
/// Devuelve el puntero a Part con los datos del PART.
/// Returns a pointer to Part with the data of the PART.
//==============================================================================
JBinaryData* JPartDataBi4::GetData()const{
  if(!Data)Run_Exceptioon("The data object is not available.");
  return(Data);
}

//==============================================================================
/// Devuelve el puntero a Part con los datos del PART.
/// Returns a pointer to Part with the data of the PART.
//==============================================================================
JBinaryData* JPartDataBi4::GetPart()const{
  if(!Part)Run_Exceptioon("PART data is not available.");
  return(Part);
}

//==============================================================================
/// Devuelve el numero de arrays en los datos del PART.
/// Returns number of arrays in PART data.
//==============================================================================
unsigned JPartDataBi4::ArraysCount()const{
  return(GetPart()->GetArraysCount());
}

//==============================================================================
/// Devuelve el nombre del array indicado.
/// Returns name of requested array.
//==============================================================================
std::string JPartDataBi4::ArrayName(unsigned num)const{
  if(num>=ArraysCount())Run_Exceptioon("Array is not available.");
  return(GetPart()->GetArray(num)->GetName());
}

//==============================================================================
/// Devuelve true cuando el array es triple.
/// Returns true when the array is triple.
//==============================================================================
bool JPartDataBi4::ArrayTriple(unsigned num)const{
  if(num>=ArraysCount())Run_Exceptioon("Array is not available.");
  return(JBinaryDataDef::TypeIsTriple(GetPart()->GetArray(num)->GetType()));
}

//==============================================================================
/// Devuelve true si existe el array indicado.
/// Returns true if the indicated array exists.
//==============================================================================
bool JPartDataBi4::ArrayExists(std::string name)const{
  return(GetPart()->GetArray(name)!=NULL);
}

//==============================================================================
/// Devuelve el puntero a Part con los datos del PART.
/// Returns a pointer to Part with the data of the PART.
//==============================================================================
JBinaryDataArray* JPartDataBi4::GetArray(std::string name)const{
  JBinaryDataArray* ar=GetPart()->GetArray(name);
  if(!ar)Run_Exceptioon(fun::PrintStr("Array \'%s\' is not available.",name.c_str()));
  return(ar);
}

//==============================================================================
/// Devuelve el puntero a Part con los datos del PART y comprueba el tipo.
/// Returns a pointer to Part with the data of the PART and checks the type.
//==============================================================================
JBinaryDataArray* JPartDataBi4::GetArray(std::string name,JBinaryDataDef::TpData type)const{
  JBinaryDataArray* ar=GetArray(name);
  if(ar->GetType()!=type)Run_Exceptioon(fun::PrintStr("Type of array \'%s\' is not %s.",name.c_str(),JBinaryDataDef::TypeToStr(type).c_str()));
  return(ar);
}

//==============================================================================
/// Devuelve el valor de Y de datos 2D.
/// Returns Y value in 2-D data.
//==============================================================================
double JPartDataBi4::Get_Particles2dPosY()const{
  double posy=0;
  if(Get_Data2d()){
    posy=DBL_MAX;
    unsigned np=Get_Npok();
    if(!np)Run_Exceptioon("Number of particles is invalid to calculates Y in 2D simulations.");
    if(Get_PosSimple()){
      tfloat3 *pos=new tfloat3[np];
      Get_Pos(np,pos);
      posy=pos[0].y;
      delete[] pos;
    }
    else{
      tdouble3 *posd=new tdouble3[np];
      Get_Posd(np,posd);
      posy=posd[0].y;
      delete[] posd;
    }
  }
  return(posy);
}