INT: Integration of functionality produced by The Environmental Hydraulics Institute "IHCantabria" (http://www.ihcantabria.com/en/)

Capabilities include:
- Wave generation
  - Solitary wave using Boussinesq theory
  - Cnoidal wave theory
  - StokesI, StokesII, StokesV wave theory

- Active wave absorption at the inflow/outflow boundaries based on
  shallow water theory

Authors:
- Javier Lopez Lara (jav.lopez@unican.es)
- Gabriel Barajas (barajasg@unican.es)
- Inigo Losada (losadai@unican.es)

integration/genAbs/allMake

@@ -0,0 +1,33 @@
+#!/bin/bash
+
+if [ $WM_PROJECT == "foam" ]; then
+    ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'`
+else
+    ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'`
+fi
+
+#IHC_dbg
+echo $ofversion
+#----
+
+export OF_VERSION=$ofversion
+
+wclean all > /dev/null
+
+wmake libso waveGeneration
+
+if (( $? )) ; then 
+    echo "\n\nWave generation boundary conditions (V2.0-ESI) compilation failed"
+    exit; else
+    echo -e "\n\nIH Wave generationV2.0 boundary conditions (V2.0-ESI) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n";  
+fi
+
+wmake libso waveAbsorption
+
+if (( $? )) ; then 
+    echo "\n\nWave absorption boundary conditions (V2.0-ESI) compilation failed"
+    exit; else
+    echo -e "\n\nIH Wave absorption boundary conditions (V2.0-ESI) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n";  
+fi
+
+wclean all > /dev/null

integration/genAbs/common/calculateWaterLevel/calculatedLevelRegular.H

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if ( waveTheory_ == "StokesI" )
+    {
+        forAll(calculatedLevel, itS1)
+        {
+            calculatedLevel[itS1] = RealwaterDepth_ 
+		+ timeMult * StokesIFun :: eta 
+		(
+    		    waveHeight_, 
+	            waveKx, 
+	            xGroup[itS1], 
+                    waveKy, 
+                    yGroup[itS1], 
+                    waveOmega, 
+                    currTime, 
+                    wavePhase_
+                );
+        }
+    }
+    else if ( waveTheory_ == "StokesII" )
+    {
+        forAll(calculatedLevel, itS2)
+        {
+            calculatedLevel[itS2] = RealwaterDepth_ 
+	        + timeMult * StokesIIFun :: eta 
+		(
+                    waveHeight_, 
+                    RealwaterDepth_,  
+                    waveKx, 
+                    xGroup[itS2], 
+                    waveKy, 
+                    yGroup[itS2], 
+                    waveOmega, 
+                    currTime, 
+                    wavePhase_
+                );
+        }
+    }
+    else  if ( waveTheory_ == "StokesV" )
+    {
+        forAll(calculatedLevel, it1)
+        {
+            calculatedLevel[it1] = RealwaterDepth_ 
+	        + timeMult * stokesVFun :: eta 
+		(
+                    RealwaterDepth_,  
+		    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+                    xGroup[it1], 
+                    yGroup[it1], 
+                    currTime, 
+                    wavePhase_
+                );
+        }
+    }
+    else if ( waveTheory_ == "Cnoidal" )
+    {
+        forAll(calculatedLevel, it3)
+        {
+            calculatedLevel[it3] = RealwaterDepth_ 
+	        + timeMult * cnoidalFun :: eta 
+	        (
+		    waveHeight_, 
+                    mCnoidal_, 
+                    waveKx, 
+                    waveKy, 
+                    wavePeriod_, 
+                    xGroup[it3], 
+                    yGroup[it3], 
+                    currTime
+                );
+        }
+    }
+    else 
+    {
+        FatalError << "Wave theory not supported, use:\n" 
+	    << "StokesI, StokesII, StokesV, Cnoidal, SolitaryBoussinesq" 
+	    << exit(FatalError);
+    }

integration/genAbs/common/calculateWaterLevel/calculatedLevelRegularNormal.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if ( waveTheory_ == "StokesI" )
+    {
+            calculatedLevel = RealwaterDepth_ 
+	        + timeMult * StokesIFun :: eta 
+		(
+		    waveHeight_, 
+                    waveKx, 
+                    xGroup[1],  
+                    waveKy, 
+                    yGroup[1],  
+                    waveOmega, 
+                    currTime, 
+                    wavePhase_
+                );
+    }
+    else if ( waveTheory_ == "StokesII" )
+    {
+            calculatedLevel = RealwaterDepth_ 
+	        + timeMult * StokesIIFun :: eta 
+	        (
+                    waveHeight_, 
+                    RealwaterDepth_,  
+                    waveKx, 
+                    xGroup[1],  
+                    waveKy, 
+                    yGroup[1],  
+                    waveOmega, 
+                    currTime, 
+                    wavePhase_
+                );
+    }
+    else  if ( waveTheory_ == "StokesV" )
+    {
+            calculatedLevel = RealwaterDepth_ 
+	        + timeMult * stokesVFun :: eta 
+		(
+                    RealwaterDepth_,  
+		    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+                    xGroup[1], 
+                    yGroup[1], 
+                    currTime, 
+                    wavePhase_
+                );
+    }
+    else if ( waveTheory_ == "Cnoidal" )
+    {
+            calculatedLevel = RealwaterDepth_ 
+	        + timeMult * cnoidalFun :: eta 
+	        (
+		    waveHeight_, 
+                    mCnoidal_, 
+                    waveKx, 
+                    waveKy, 
+                    wavePeriod_, 
+                    xGroup[1], 
+                    yGroup[1], 
+                    currTime
+                );
+    }
+    else 
+    {
+        FatalError << "Wave theory not supported, use:\n" 
+	    << "StokesI, StokesII, StokesV, Cnoidal, SolitaryBoussinesq." 
+	    << exit(FatalError);
+    }

integration/genAbs/common/calculateWaterLevel/calculatedLevelSolitary.H

@@ -0,0 +1,30 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/ 
+
+   if ( waveTheory_ == "Boussinesq" )
+   {
+       forAll(calculatedLevel, it2)
+       {
+           calculatedLevel[it2] = RealwaterDepth_
+               + BoussinesqFun :: eta
+               (
+                   waveHeight_, 
+                   RealwaterDepth_,  
+                   xGroup[it2], 
+                   yGroup[it2], 
+                   waveAngle, 
+                   currTime, 
+                   X0
+               );
+        }
+   }
+   else 
+   {
+        FatalError << "Wave theory not supported, use:\n" 
+	    << "Boussinesq" << exit(FatalError);
+   }

integration/genAbs/common/checks/checkInputErrorsRegular.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+if (nPaddles_ < 1)
+{
+    FatalError << "Check nPaddles value." << exit(FatalError);
+}
+
+if ( nPaddles_ > 1 )
+{
+    nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan );
+    reduce(nPaddles_, minOp<label>());
+}
+
+if ( waveHeight_ < 0.0 )
+{
+    FatalError << "Check wave height value." << exit(FatalError);
+}
+
+if ( waveTheory_ == "aaa" ) 
+{
+    FatalError << "Wave theory not specified." << exit(FatalError);
+}
+else if ( waveTheory_ == "StokesI" || waveTheory_ == "StokesII" )
+{
+    if ( wavePeriod_ <= 0.0 )
+    {
+        FatalError << "Check wave period value." << exit(FatalError);
+    }
+
+    if ( waveLength_ <= 0.0 )
+    {
+        word generation = "a";
+
+        waveLength_ = StokesIFun::waveLength (RealwaterDepth_, wavePeriod_);
+
+        waveK = 2.0*PI()/waveLength_;
+
+        if ( waveK*RealwaterDepth_ > PI() )
+        {
+            generation = "Deep";
+        }
+        else if ( waveK*RealwaterDepth_ < PI()/10.0 )
+        {
+            generation = "Shallow";
+        }
+        else
+        {
+            generation = "Intermediate";
+        }
+
+        Info << "\nIH Wave Generation BC" << endl;
+        Info << "	Wave theory: " << waveTheory_ << endl;
+        Info << "	H: " << waveHeight_ << endl;
+        Info << "	T: " << wavePeriod_ << endl;
+        Info << "	h: " << RealwaterDepth_ << endl;
+        Info << "	waveLength: " << waveLength_ << endl;
+        Info << "	Direction: " << waveDir_ << "º" << endl;
+        Info << "	Generation in: " << generation << " waters." << endl;
+        Info << "	Relative depth (kh): " << waveK*RealwaterDepth_ 
+	     << "\n\n" << endl;
+    }
+}
+else if ( waveTheory_ == "StokesV" )
+{
+    if ( wavePeriod_ <= 0.0 )
+    {
+        FatalError << "Check wave period value." << exit(FatalError);
+    }
+
+    if ( lambdaStokesV_ <= 0.0 )
+    {
+        word generation = "a";
+
+        scalar f1;
+        scalar f2;
+
+        stokesVFun :: StokesVNR ( waveHeight_, 
+            RealwaterDepth_,  
+            wavePeriod_, 
+            &waveK, 
+            &lambdaStokesV_, 
+            &f1, 
+            &f2 
+            );
+
+        Info << "f1 residual " << f1 << endl;
+        Info << "f2 residual " << f2 << endl;
+
+        if ( f1 > 0.001 || f2 > 0.001 )
+        {
+	    FatalError << "No convergence for Stokes V wave theory.\n" 
+	        << exit(FatalError);
+        }
+
+        waveLength_ = 2.0*PI()/waveK;
+
+        if ( waveK*RealwaterDepth_ > PI() )
+        {
+            generation = "Deep";
+        }
+        else if ( waveK*RealwaterDepth_ < PI()/10.0 )
+        {
+            generation = "Shallow";
+        }
+        else
+        {
+            generation = "Intermediate";
+        }
+
+        Info << "\nIH Wave Generation BC" << endl;
+        Info << "	Wave theory: " << waveTheory_ << endl;
+        Info << "	H: " << waveHeight_ << endl;
+        Info << "	T: " << wavePeriod_ << endl;
+        Info << "	h: " << RealwaterDepth_ << endl;
+        Info << "	L: " << waveLength_ << endl;
+        Info << "	Lambda: " << lambdaStokesV_ << endl;
+        Info << "	Direction: " << waveDir_ << "º"  << endl;
+        Info << "	Generation in: " << generation << " waters." << endl;
+        Info << "	Relative depth (kh): " << waveK*RealwaterDepth_ 
+	     << "\n\n" << endl;
+    }
+}
+else if ( waveTheory_ == "Cnoidal" )
+{
+    if ( wavePeriod_ <= 0.0 )
+    {
+        FatalError << "Check wave period value." << exit(FatalError);
+    }
+
+    if ( mCnoidal_ <= 0.0 )
+    {
+        cnoidalFun :: calculations ( waveHeight_, 
+            RealwaterDepth_,  
+            wavePeriod_, 
+            &mCnoidal_, 
+            &waveLength_);
+
+        Info << "\nIH Wave Generation BC" << endl;
+        Info << "	Wave theory: " << waveTheory_ << endl;
+        Info << "	H: " << waveHeight_ << endl;
+        Info << "	T: " << wavePeriod_ << endl;
+        Info << "	h: " << RealwaterDepth_ << endl;
+        Info << "	L: " << waveLength_ << endl;
+        Info << "	m parameter: " << mCnoidal_ << endl;
+        Info << "	Direction: " << waveDir_ << "º" <<  "\n\n" << endl;
+    }
+}
+else
+{
+    FatalError << "Wave theory not supported, use:\n" 
+        << "StokesI, StokesII, StokesV, Cnoidal." << exit(FatalError);
+}
+

integration/genAbs/common/checks/checkInputErrorsSolitary.H

@@ -0,0 +1,41 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (nPaddles_ < 1)
+    {
+        FatalError << "Check nPaddles value." << exit(FatalError);
+    }
+
+    if ( nPaddles_ > 1 )
+    {
+        nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan );
+        reduce(nPaddles_, minOp<label>());
+    }
+
+    if ( waveHeight_ < 0.0 )
+    {
+        FatalError << "Check wave height value." << exit(FatalError);
+    }
+
+    if ( waveTheory_ == "aaa" )
+    {
+        FatalError << "Wave theory not specified." << exit(FatalError);
+    }
+    else if ( waveTheory_ == "Boussinesq" )
+    {
+        Info << "\nIH Wave Generation BC" << endl;
+        Info << "Wave theory: " << waveTheory_ << endl;
+        Info << "H: " << waveHeight_ << endl;
+        Info << "h: " << RealwaterDepth_ << endl;
+        Info << "Direction: " << waveDir_ << "º" <<  "\n\n" << endl;
+    }
+    else
+    {
+        FatalError << "Wave theory not supported, use:\n" 
+	    << "Boussinesq" << exit(FatalError);
+    }

integration/genAbs/common/memberFun.H

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    //- Return Pi
+    scalar PI()
+    {
+        #if OFVERSION >= 200
+            const scalar PI = constant::mathematical::pi;
+        #else
+            const scalar PI = mathematicalConstant::pi;
+        #endif
+        return PI;
+    }
+
+    //- Return alphaName
+    word alphaName()
+    {
+        #if OFVERSION >= 230
+            const word an = "alpha.water";
+        #else
+            const word an = "alpha1";
+        #endif
+        return an;
+    }
+
+    //- Return prevalent direction - X or Y to sort groups of columns
+    scalarField patchDirection ( vector span, scalar* dMin, scalar* dSpan )
+    {
+        const pointField& points = patch().patch().localPoints();
+        if (span[0] >= span[1]) 
+        {
+            *dMin = gMin(points.component(0));
+            *dSpan = gMax(points.component(0)) - *dMin;
+            return patch().Cf().component(0); 
+        }
+        else
+        {
+            *dMin = gMin(points.component(1));
+            *dSpan = gMax(points.component(1)) - *dMin;
+            return patch().Cf().component(1);
+        }
+    }
+
+    //- Return Zmax and Zmin of the patch faces
+    void faceBoundsZ(scalarField* zSup, scalarField* zInf)
+    {
+        const label nF = patch().faceCells().size();
+        scalarField zMax = Foam::scalarField(nF, -9999.0);
+        scalarField zMin = Foam::scalarField(nF, 9999.0);
+            
+        const faceList& faces = this->patch().patch().localFaces();
+        const List<vector>& points = this->patch().patch().localPoints();
+
+        forAll( faces, faceI )
+        {
+            const face& f = faces[faceI];
+            forAll(f, point)
+            {
+                scalar auxiliar = points[f[point]].component(2);
+
+                zMax[faceI] = max(zMax[faceI], auxiliar);
+                zMin[faceI] = min(zMin[faceI], auxiliar);
+            }
+        }
+
+        *zSup = zMax;
+        *zInf = zMin;
+    }
+
+    //- Calculate water levels on each paddle (same zSpan)
+    scalarList calcWL 
+        ( 
+	    scalarField alphaCell, 
+	    labelList cellGroup, 
+	    scalar zSpan 
+	)
+    {
+        scalarList groupTotalArea (nPaddles_,0.0); 
+        scalarList groupWaterArea (nPaddles_,0.0); 
+        scalarList heights (nPaddles_,0.0);
+
+        const scalarField faceSurface = patch().magSf();
+
+        forAll( patch().faceCells(), index ) 
+        {
+            groupTotalArea[cellGroup[index]-1] 
+	        += faceSurface[index];
+            groupWaterArea[cellGroup[index]-1] 
+	        += faceSurface[index]*alphaCell[index];
+        } 
+
+        for (label i=0; i<=gMax(cellGroup)-1; i++)
+        {
+            reduce(groupTotalArea[i], sumOp<scalar>());
+            reduce(groupWaterArea[i], sumOp<scalar>());
+            heights[i] = groupWaterArea[i]/groupTotalArea[i]*zSpan;
+        }
+
+        return heights;
+    }
+
+    //- Calculate mean velocities on each paddle
+    void calcUV 
+        ( 
+	    scalarField alphaCell, 
+	    labelList cellGroup, 
+	    scalarField UxCell, 
+            scalarList* Ux, 
+	    scalarField UyCell, 
+	    scalarList* Uy 
+	)
+    {
+       scalarList groupWaterArea (nPaddles_,0.0);
+       scalarList groupUx (nPaddles_,0.0);
+       scalarList groupUy (nPaddles_,0.0);
+
+       const scalarField faceSurface = patch().magSf();
+
+       forAll( patch().faceCells(), index ) 
+       {
+           groupWaterArea[cellGroup[index]-1] += 
+               faceSurface[index]*alphaCell[index];
+
+           groupUx[cellGroup[index]-1] += 
+               UxCell[index]*alphaCell[index]*faceSurface[index];
+           groupUy[cellGroup[index]-1] += 
+               UyCell[index]*alphaCell[index]*faceSurface[index];
+       } 
+
+       for (label i=0; i<=nPaddles_-1; i++)
+       {
+           reduce(groupWaterArea[i], sumOp<scalar>());
+            reduce(groupUx[i], sumOp<scalar>());
+            reduce(groupUy[i], sumOp<scalar>());
+
+            groupUx[i] = groupUx[i]/groupWaterArea[i];
+            groupUy[i] = groupUy[i]/groupWaterArea[i];
+        }
+
+        *Ux = groupUx;
+        *Uy = groupUy;
+    }
+
+    //- Mean of a scalarList
+    scalar meanSL ( scalarList lst )
+    {
+        scalar aux = 0.0;
+            
+        forAll(lst,i)
+        {
+            aux += lst[i];
+        }
+            
+        return aux/scalar(lst.size());
+    }
+
+    //- Reduce number of paddles if too high (paddles without cells)
+    label decreaseNPaddles 
+        ( 
+	    label np, 
+	    scalarField patchD, 
+	    scalar dMin, 
+	    scalar dSpan 
+	)
+    {
+        scalarList dBreakPoints(np+1, 0.0); 
+
+        while(1)
+        {
+            scalarList usedGroup (np,0.0); 
+
+            for (label i=0; i<=np; i++)
+            {
+                dBreakPoints[i] = dMin + dSpan/(np)*i; 
+            }
+
+            forAll(patch().faceCells(), patchCells) 
+            {
+                for (label j=0; j<=np; j++) 
+                {
+                    if ( (patchD[patchCells]>=dBreakPoints[j]) && 
+                        (patchD[patchCells]<dBreakPoints[j+1]) )
+                    {
+                        usedGroup[j] = 1.0;
+                        continue;
+                    }
+                }      
+            }
+
+            reduce(usedGroup, maxOp<scalarList>());
+
+            if ( gMin(usedGroup) < 0.5 )
+            {
+                Info
+                    << "Reduced nPaddles from " << np
+                    << " to " << np-1 << endl;
+                np -= 1;
+            }
+            else
+            {
+                break;
+            }
+        }
+        return np;
+    }
+
+    //- Own convenient definition of atan
+    scalar arcTan ( scalar x, scalar y )
+    {
+        scalar A = 0.0;
+
+        if( x > 0.0 )
+        {
+            A = atan(y/x);
+        }
+        else if ( y >= 0.0 && x < 0.0 )
+        {
+            A = PI() + atan(y/x);
+        }
+        else if ( y < 0.0 && x < 0.0 )
+        {
+            A = -PI() + atan(y/x);
+        }
+        else if ( y > 0.0 && x == 0.0 )
+        {
+            A = PI()/2.0;
+        }        
+        else if ( y < 0.0 && x == 0.0 )
+        {
+            A = -PI()/2.0;
+        }
+        else if ( y == 0.0 && x == 0.0 )
+        {
+            A = 0;
+        }
+
+        return A;
+    }
+
+    //- Calculate mean horizontal angle for each paddle
+    scalarList meanPatchDirs ( labelList cellGroup )
+    {
+        const vectorField nVecCell = patch().nf();
+
+        scalarList numAngles (nPaddles_, 0.0);
+        scalarList meanAngle (nPaddles_, 0.0);
+
+        forAll(patch().faceCells(), patchCells) 
+        {
+            meanAngle[cellGroup[patchCells]-1]
+                += arcTan( nVecCell[patchCells].component(0), 
+                nVecCell[patchCells].component(1));
+            numAngles[cellGroup[patchCells]-1] += 1.0;
+        }
+
+        for (label i=0; i<=nPaddles_-1; i++)
+        {
+            reduce(meanAngle[i], sumOp<scalar>());
+            reduce(numAngles[i], sumOp<scalar>());
+            meanAngle[i] = meanAngle[i]/numAngles[i] + PI();
+        }
+        return meanAngle;
+    }
+
+    //- Calculate z-bounds for each paddle
+    scalarList zSpanList 
+        ( 
+	    labelList cellGroup, 
+	    scalarField zInf, 
+	    scalarField zSup 
+	)
+    {
+        scalarList zMaxL (nPaddles_, -9999.0);
+        scalarList zMinL (nPaddles_,  9999.0);
+
+        forAll(patch().faceCells(), patchCells) 
+        {
+            zMaxL[cellGroup[patchCells]-1] = 
+	        max(zMaxL[cellGroup[patchCells]-1], zSup[patchCells] );
+            zMinL[cellGroup[patchCells]-1] = 
+                min(zMinL[cellGroup[patchCells]-1], zInf[patchCells] );
+        }
+
+        for (label i=0; i<=nPaddles_-1; i++)
+        {
+            reduce(zMaxL[i], maxOp<scalar>());
+            reduce(zMinL[i], minOp<scalar>());
+        }
+        return zMaxL-zMinL;
+    }
+
+    //- In-line print for scalarLists
+    label inlinePrint ( std::string name, scalarList SL )
+    {
+        Info << name << " " << SL.size() << "( ";
+        forAll(SL, i)
+        {
+            Info << SL[i] << " ";
+        }
+        Info << ")\n" << endl;
+        return 0;
+    }
+
+    //- Simple linear interpolation
+    scalar interpolation 
+        (
+	    scalar x1, 
+	    scalar x2, 
+            scalar y1, 
+            scalar y2, 
+            scalar xInt 
+        )
+    {
+        scalar yInt = y1 + (y2-y1)/(x2-x1)*(xInt-x1);
+        return yInt;
+    }
+
+    //- Limit angle between -pi/2 and pi/2
+    scalar limAngle (scalar ang)
+    {
+        ang = abs(ang);
+
+        while (ang >= 2.0*PI())
+        {
+            ang -= 2.0*PI();
+        }
+
+        if ( ang >= PI()/2.0 && ang <= 3.0*PI()/2.0 )
+        {
+            return PI()/2.0;
+        }
+        else
+        {
+            return ang;
+        }
+    }
+
+    //- Interpolation BO style
+    scalar interpolatte 
+        ( 
+            scalar x1, 
+            scalar x2, 
+            scalar y1, 
+            scalar y2, 
+            scalar xInt 
+        )
+    {
+        scalar yInt = (y1 + y2)*0.5;
+        return yInt;
+    }

integration/genAbs/common/waveFun.C

@@ -0,0 +1,980 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#include "waveFun.H"
+#include <math.h>
+#include <stdlib.h> 
+#include <stdio.h>
+
+namespace otherFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define grav 9.81
+
+    double interpolation (double x1, double x2, double y1, double y2, double xInt)
+    {
+
+        double yInt = y1 + (y2-y1)/(x2-x1)*(xInt-x1);
+
+        return yInt;
+    }
+}
+
+namespace StokesIFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define G 9.81
+
+    double waveLength (double h, double T)
+    {
+        double L0 = G*T*T/(2.0*PII);
+        double L = L0;
+
+        for(int i=1; i<=100; i++)
+        {
+            L = L0*tanh(2.0*PII*h/L);
+        }
+
+        return L;
+    }
+    
+    double eta (double H, double Kx, double x, double Ky, double y, double omega, double t, double phase)
+    {
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+        
+        double sup = H*0.5*cos(faseTot);
+        
+        return sup;
+    }
+
+    double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z)
+    {
+        double k = sqrt(Kx*Kx + Ky*Ky);
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+
+        double velocity = H*0.5*omega*cos(faseTot)*cosh(k*z)/sinh(k*h);
+
+        return velocity;
+    }
+    
+    double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z)
+    {
+        double k = sqrt(Kx*Kx + Ky*Ky);
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+
+        double velocity = H*0.5*omega*sin(faseTot)*sinh(k*z)/sinh(k*h);
+
+        return velocity;
+    }   
+}
+
+namespace StokesIIFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define G 9.81 
+    
+    double eta (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase)
+    {
+        double k = sqrt(Kx*Kx + Ky*Ky);
+        double sigma = tanh(k*h);
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+        
+        double sup = H*0.5*cos(faseTot) + k*H*H/4.0*(3.0-sigma*sigma)/(4.0*pow(sigma,3))*cos(2.0*faseTot);
+        
+        return sup;
+    }
+    
+    double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z)
+    {
+        double k = sqrt(Kx*Kx + Ky*Ky);
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+
+        double velocity = H*0.5*omega*cos(faseTot)*cosh(k*z)/sinh(k*h) + 3.0/4.0*H*H/4.0*omega*k*cosh(2.0*k*z)/pow(sinh(k*h),4)*cos(2.0*faseTot);
+
+        return velocity;
+    }
+    
+    double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z)
+    {
+        double k = sqrt(Kx*Kx + Ky*Ky);
+        double faseTot = Kx*x + Ky*y - omega*t + phase;
+
+        double velocity = H*0.5*omega*sin(faseTot)*sinh(k*z)/sinh(k*h) + 3.0/4.0*H*H/4.0*omega*k*sinh(2.0*k*z)/pow(sinh(k*h),4)*sin(2.0*faseTot);
+
+        return velocity;
+    }      
+}
+
+namespace Elliptic
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define ITER 25
+
+    void ellipticIntegralsKE (double m, double* K, double* E)
+    {
+        long double a,aOld,g,gOld,aux,sum;
+
+        if ( m == 0.0 ) 
+        {
+            *K = PII/2.;
+            *E = PII/2.;
+            return;
+        }
+
+        a = 1.0L;
+        g = sqrt(1.0L - m);
+        aux = 1.0L;
+        sum = 2.0L - m;
+
+        while (1) 
+        {
+            gOld = g;
+            aOld = a;
+            a = 0.5L * (gOld + aOld);
+            g = gOld * aOld;
+            aux += aux;
+            sum -= aux * (a * a - g);
+
+            if ( fabs(aOld - gOld) <= (aOld * 1.e-22) )
+            { 
+                break;
+            }
+
+            g = sqrt(g);
+        }
+
+        *K = (double) (PII/2. / a);
+        *E = (double) (PII/4. / a * sum);
+        return;
+    }
+
+    double JacobiAmp (double u, double m)
+    {
+        long double a[ITER+1], g[ITER+1], c[ITER+1];
+        long double aux, amp;
+        int n;
+
+        m = fabsl(m);
+
+        if ( m == 0.0 ) 
+        {
+            return u;
+        }
+        
+        if ( m == 1. ) 
+        {
+            return 2. * atan( exp(u) ) - PII/2.;
+        }
+
+        a[0] = 1.0L;
+        g[0] = sqrtl(1.0L - m);
+        c[0] = sqrtl(m);
+
+        aux = 1.0L;
+
+        for (n = 0; n < ITER; n++) 
+        {
+            if ( fabsl(a[n] - g[n]) < (a[n] * 1.e-22) ) 
+            {
+                break;
+            }
+
+            aux += aux;
+            a[n+1] = 0.5L * (a[n] + g[n]);
+            g[n+1] = sqrtl(a[n] * g[n]);
+            c[n+1] = 0.5L * (a[n] - g[n]);
+        }
+
+        amp = aux * a[n] * u;
+
+        for (; n > 0; n--) 
+        {
+            amp = 0.5L * ( amp + asinl( c[n] * sinl(amp) / a[n]) );
+        }
+
+        return (double) amp; 
+    }
+
+    void JacobiSnCnDn (double u, double m, double* Sn, double* Cn, double* Dn)
+    {
+        double amp = Elliptic::JacobiAmp( u, m );
+
+        *Sn = sin( amp );
+        *Cn = cos( amp );
+        *Dn = sqrt(1.0 - m * sin( amp ) * sin( amp ));
+        return;
+    }
+
+}
+
+namespace cnoidalFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define G 9.81
+
+    double eta (double H, double m, double kx, double ky, double T, double x, double y, double t)
+    {
+        double K, E;
+        Elliptic::ellipticIntegralsKE(m, &K, &E);
+
+        double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T);
+
+        double sn, cn, dn;
+        Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn);
+
+        double etaCnoidal = H*((1.0-E/K)/m - 1.0 + pow(cn,2));
+
+        return etaCnoidal;
+    }
+
+    double etaCnoidal1D (double H, double m, double t, double T)
+    {
+        double K, E;
+        Elliptic::ellipticIntegralsKE(m, &K, &E);
+
+        double uCnoidal = -2.0*K*(t/T);
+
+        double sn, cn, dn;
+        Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn);
+
+        double etaCnoidal = H*((1.0-E/K)/m - 1.0 + pow(cn,2));
+
+        return etaCnoidal;
+    }
+
+    double etaMeanSq (double H, double m, double T)
+    {
+        double eta = 0;
+        double etaSumSq = 0;
+
+        for(int i=0; i<1000; i++)
+        {
+            eta = etaCnoidal1D(H, m, i*T/(1000.0), T);
+            etaSumSq += eta*eta;
+        }
+
+        etaSumSq /= 1000.0;
+        return etaSumSq;
+    }
+
+    double d1EtaDx (double H, double m, double uCnoidal, double L, double K, double E)
+    {
+        double dudx = 0;
+        dudx = 2.0*K/L;
+
+        double sn, cn, dn;
+        Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn);
+
+	double deriv = -2.0*H*cn*dn*sn*dudx;
+
+        return deriv;
+    }
+
+    double d2EtaDx (double H, double m, double uCnoidal, double L, double K, double E)
+    {
+        double dudxx = 0;
+        dudxx = 4.0*K*K/L/L;
+
+        double sn, cn, dn;
+        Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn);
+
+        double deriv = 2.0*H*(dn*dn*sn*sn - cn*cn*dn*dn + m*cn*cn*sn*sn)*dudxx;
+
+        return deriv;
+    }
+
+    double d3EtaDx (double H, double m, double uCnoidal, double L, double K, double E)
+    {
+        double dudxxx = 0;
+        dudxxx = 8.0*K*K*K/L/L/L;
+
+        double sn, cn, dn;
+        Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn);
+
+        double deriv = 8.0*H*( cn*sn*dn*dn*dn*(-4.0 -2.0*m) + 4.0*m*cn*sn*sn*sn*dn -2.0*m*cn*cn*cn*sn*dn )*dudxxx;
+
+        return deriv;
+    }
+
+    int calculations (double H, double d, double T, double* mOut, double* LOut)
+    {
+        double mTolerance = 0.0001;
+        double mElliptic = 0.5;
+        double LElliptic = 0;
+        double phaseSpeed = 0;
+
+        double mError = 0.0;
+        double mMinError = 999;
+
+        while (mElliptic < 1.0)
+        {
+            double KElliptic, EElliptic;
+            Elliptic::ellipticIntegralsKE(mElliptic, &KElliptic, &EElliptic);
+
+            LElliptic = KElliptic*sqrt(16.0*pow(d,3)*mElliptic/(3.0*H));        
+
+            phaseSpeed = sqrt(G*d)*(1.0 - H/d/2.0 + H/d/mElliptic*(1.0-3.0/2.0*EElliptic/KElliptic));
+
+            mError = fabs(T-LElliptic/phaseSpeed);
+
+            if (mError <= mMinError)
+            {
+                *mOut = mElliptic;
+                *LOut = LElliptic;
+                mMinError = mError;
+            }
+
+            mElliptic += mTolerance;
+        }
+
+        return 0;
+    }
+
+    double U (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z)
+    {
+        double K, E;
+        Elliptic::ellipticIntegralsKE(m, &K, &E);
+
+        double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T);
+        double k = sqrt(kx*kx + ky*ky);
+        double L = 2.0*PII/k;
+        double c = L/T;
+
+        double etaCN = eta(H, m, kx, ky, T, x, y, t);
+        double etaXX = d2EtaDx(H, m, uCnoidal, L, K, E);
+        double etaMS = etaMeanSq(H, m, T);
+
+        double velocity = c*etaCN/h - c*(etaCN*etaCN/h/h + etaMS*etaMS/h/h) + 1.0/2.0*c*h*(1.0/3.0 - z*z/h/h)*etaXX;
+
+        return velocity;
+    }
+
+    double W (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z)
+    {
+        double K, E;
+        Elliptic::ellipticIntegralsKE(m, &K, &E);
+
+        double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T);
+        double k = sqrt(kx*kx + ky*ky);
+        double L = 2.0*PII/k;
+        double c = L/T;
+
+        double etaCN = eta(H, m, kx, ky, T, x, y, t);
+        double etaX = d1EtaDx(H, m, uCnoidal, L, K, E);
+        double etaXXX = d3EtaDx(H, m, uCnoidal, L, K, E);
+
+
+        double velocity = -c*z*( etaX/h*(1.0-2.0*etaCN/h) + 1.0/6.0*h*(1.0-z*z/h/h)*etaXXX ) ;
+
+        return velocity;
+    }
+}
+
+namespace stokesVFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define G 9.81
+
+    double A11 (double h, double k)
+    {
+        double s = sinh(k*h);
+
+        double A = 1.0/s;
+
+        return A;
+    }
+
+    double A13 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            -pow(c,2)*(5.0*pow(c,2)+1.0)
+                /
+            (8.0*pow(s,5));
+
+        return A;
+    }
+ 
+    double A15 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            -(1184.0*pow(c,10)-1440.0*pow(c,8)-1992.0*pow(c,6)+2641.0*pow(c,4)-249.0*pow(c,2)+18)
+                /
+            (1536.0*pow(s,11));
+
+        return A;
+    }   
+
+    double A22 (double h, double k)
+    {
+        double s = sinh(k*h);
+
+        double A = 
+            3.0
+                /
+            (8.0*pow(s,4));
+
+        return A;
+    }
+
+    double A24 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            (192.0*pow(c,8)-424.0*pow(c,6)-312.0*pow(c,4)+480.0*pow(c,2)-17)
+                /
+            (768.0*pow(s,10));
+
+        return A;
+    }
+
+    double A33 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            (13.0-4.0*pow(c,2))
+                /
+            (64.0*pow(s,7));
+
+        return A;
+    }
+
+    double A35 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            (512.0*pow(c,12)+4224.0*pow(c,10)-6800.0*pow(c,8)-12808.0*pow(c,6)+16704.0*pow(c,4)-3154.0*pow(c,2)+107.0)
+                /
+            (4096.0*pow(s,13)*(6.0*pow(c,2)-1.0));
+
+        return A;
+    }
+
+    double A44 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            (80.0*pow(c,6)-816.0*pow(c,4)+1338.0*pow(c,2)-197.0)
+                /
+            (1536.0*pow(s,10)*(6.0*pow(c,2)-1.0));
+
+        return A;
+    }
+
+    double A55 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double A = 
+            -(2880.0*pow(c,10)-72480.0*pow(c,8)+324000.0*pow(c,6)-432000.0*pow(c,4)+163470.0*pow(c,2)-16245.0)
+                /
+            (61440.0*pow(s,11)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0));
+
+        return A;
+    }
+
+    double B22 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (2.0*pow(c,2)+1)*c
+                /
+            (4.0*pow(s,3));
+
+        return B;
+    }
+
+    double B24 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (272.0*pow(c,8)-504.0*pow(c,6)-192.0*pow(c,4)+322.0*pow(c,2)+21.0)*c
+                /
+            (384.0*pow(s,9));
+
+        return B;
+    }
+
+    double B33 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (8.0*pow(c,6)+1.0)*3.0
+                /
+            (64.0*pow(s,6));
+
+        return B;
+    }
+
+    double B33k (double h, double k) // d B33 / d k
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double sk = h*s;
+        double ck = h*c;
+
+        double B = 
+            9.0*pow(c,5)*ck/(4.0*pow(s,6))- 
+            (9.0*(8.0*pow(c,6)+1.0))/(32.0*pow(s,7))*sk;
+
+        return B;
+    }
+
+    double B35 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8)-21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81)
+                /
+            (12288.0*pow(s,12)*(6.0*pow(c,2)-1.0));
+
+        return B;
+    }
+
+    double B35k (double h, double k) // d B35 / d k
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double sk = h*s;
+        double ck = h*c;
+
+        double B = 
+            (14.0*88128.0*pow(c,13)*ck-12.0*208224.0*pow(c,11)*ck
+            +10.0*70848.0*pow(c,9)*ck+8.0*54000.0*pow(c,7)*ck
+            -6.0*21816.0*pow(c,5)*ck+4.0*6264.0*pow(c,3)*ck
+            -2.0*54.0*c*ck)
+            /(12288.0*pow(s,12)*(6.0*pow(c,2)-1.0))
+            -(88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8) 
+            -21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81.0)*12.0 
+            /(12288.0*pow(s,13)*(6.0*pow(c,2)-1.0))*sk 
+            -(88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8) 
+            -21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81.0)*12.0*c*ck 
+            /(12288.0*pow(s,12)*pow((6.0*pow(c,2)-1.0),2));
+
+        return B;
+    }
+
+    double B44 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (768.0*pow(c,10)-448.0*pow(c,8)-48.0*pow(c,6)+48.0*pow(c,4)+106.0*pow(c,2)-21.0)*c
+                /
+            (384.0*pow(s,9)*(6.0*pow(c,2)-1.0));
+
+        return B;
+    }
+
+    double B55 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double B = 
+            (192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10)-7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0)
+                /
+            (12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0));
+
+        return B;
+    }
+
+    double B55k (double h, double k) // d B55 / d k
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double sk = h*s;
+        double ck = h*c;
+
+        double B = 
+            (16.0*192000.0*pow(c,15)*ck-14.0*262720.0*pow(c,13)*ck
+            +12.0*83680.0*pow(c,11)*ck+10.0*20160.0*pow(c,9)*ck
+            -8.0*7280.0*pow(c,7)*ck+6.0*7160.0*pow(c,5)*ck
+            -4.0*1800.0*pow(c,3)*ck-2.0*1050.0*pow(c,1)*ck)
+            /(12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0))
+            -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) 
+            -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0)*10.0
+            /(12288.0*pow(s,11)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0))*sk
+            -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) 
+            -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0)
+            *12.0*c*ck 
+            /(12288.0*pow(s,10)*pow((6.0*pow(c,2)-1.0),2)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0))
+            -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) 
+            -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0)
+            *(32.0*pow(c,3)-22.0*c)*ck
+            /(12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*pow((8.0*pow(c,4)-11.0*pow(c,2)+3.0),2));
+
+        return B;
+    }
+
+    double C1 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double C = 
+            (8.0*pow(c,4)-8.0*pow(c,2)+9.0)
+                /
+            (8.0*pow(s,4));
+
+        return C;
+    }
+
+    double C1k (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double sk = h*s;
+        double ck = h*c;
+
+        double C = 
+            (4.0*8.0*pow(c,3)*ck-2.0*8.0*c*ck)/(8.0*pow(s,4))
+            -(8.0*pow(c,4)-8.0*pow(c,2)+9.0)*4.0*sk/(8.0*pow(s,5));
+
+        return C;
+    }
+
+    double C2 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double C = 
+            (3840.0*pow(c,12)-4096.0*pow(c,10) + 2592.0*pow(c,8)-1008.0*pow(c,6)+5944.0*pow(c,4)-1830.0*pow(c,2)+147.0) // - 2592
+                /
+            (512.0*pow(s,10)*(6.0*pow(c,2)-1.0));
+
+        return C;
+    }
+
+    double C2k (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double sk = h*s;
+        double ck = h*c;
+
+        double C = 
+            (12.0*3840.0*pow(c,11)*ck-10.0*4096.0*pow(c,9)*ck
+            +8.0*2592.0*pow(c,7)*ck-6.0*1008.0*pow(c,5)*ck+
+            4.0*5944.0*pow(c,3)*ck-2.0*1830.0*c*ck)
+            /(512.0*pow(s,10)*(6.0*pow(c,2)-1.0))
+            -(3840.0*pow(c,12)-4096.0*pow(c,10)+2592.0*pow(c,8)-1008.0*pow(c,6)+
+            5944.0*pow(c,4)-1830.0*pow(c,2)+147.0)*10.0*sk/
+            (512.0*pow(s,11)*(6.0*pow(c,2)-1.0))
+            -(3840.0*pow(c,12)-4096.0*pow(c,10)+2592.0*pow(c,8)-1008.0*pow(c,6)+
+            5944.0*pow(c,4)-1830.0*pow(c,2)+147.0)*12.0*c*ck
+            /(512.0*pow(s,10)*pow((6.0*pow(c,2)-1.0),2));
+
+        return C;
+    }
+
+    double C3 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double C = 
+            (-1.0)
+                /
+            (4.0*s*c);
+
+        return C;
+    }
+
+    double C4 (double h, double k)
+    {
+        double s = sinh(k*h);
+        double c = cosh(k*h);
+
+        double C = 
+            (12.0*pow(c,8)+36.0*pow(c,6)-162.0*pow(c,4)+141.0*pow(c,2)-27.0)
+                /
+            (192.0*c*pow(s,9));
+
+        return C;
+    }
+
+    int StokesVNR (double H, double d, double T, double* kOut, double* LambdaOut, double* f1Out, double* f2Out )
+    {
+        double f1 = 1;
+        double f2 = 1;
+
+        double k = 2.0*PII/(sqrt(G*d)*T);
+        double Lambda = H/2.0*k;
+
+        double Bmat11, Bmat12, Bmat21, Bmat22;
+        double b33, b35, b55, b33k, b35k, b55k;
+        double c1, c2, c1k, c2k;
+        double kPr, LambdaPr;
+
+        int n = 0;
+
+        while ( (fabs(f1)>1.0e-12 || fabs(f2)>1.0e-12) && n<10000 )
+        {
+            b33 = B33(d, k);
+            b35 = B35(d, k);
+            b55 = B55(d, k);
+            c1 = C1(d, k);
+            c2 = C2(d, k);
+
+            b33k = B33k(d, k);
+            b35k = B35k(d, k);
+            b55k = B55k(d, k);
+            c1k = C1k(d, k);
+            c2k = C2k(d, k);
+
+            Bmat11 = 2.0*PII/(pow(k,2)*d)*(Lambda+pow(Lambda,3)*b33+pow(Lambda,5)*(b35+b55))
+                    -2.0*PII/(k*d)*(pow(Lambda,3)*b33k+pow(Lambda,5)*(b35k+b55k));
+
+            Bmat12 = -2.0*PII/(k*d)*
+                     (1.0+3.0*pow(Lambda,2)*b33+5.0*pow(Lambda,4)*(b35+b55));
+
+            Bmat21 = -d/(2.0*PII)*tanh(k*d)*(1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2)
+                    -k*d/(2.0*PII)*(1.0-pow((tanh(k*d)),2))*d*
+                    (1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2)
+                    -k*d/(2.0*PII)*tanh(k*d)
+                    *(pow(Lambda,2)*c1k+pow(Lambda,4)*c2k);
+
+            Bmat22 = -k*d/(2.0*PII)*tanh(k*d)
+                     *(2.0*Lambda*c1+4.0*pow(Lambda,3)*c2);
+
+            f1 = PII*H/d-2.0*PII/(k*d)*(Lambda+pow(Lambda,3)*b33+pow(Lambda,5)*(b35+b55));
+
+            f2 = (2.0*PII*d)/(G*pow(T,2))-k*d/(2.0*PII)*tanh(k*d)
+                *(1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2);
+
+            LambdaPr = (f1*Bmat21-f2*Bmat11)/(Bmat11*Bmat22-Bmat12*Bmat21);
+            kPr = (f2*Bmat12-f1*Bmat22)/(Bmat11*Bmat22-Bmat12*Bmat21);
+
+            Lambda += LambdaPr;
+            k += kPr;
+
+            n++;
+        }
+
+        *kOut = k;
+        *LambdaOut = Lambda;
+
+        *f1Out = fabs(f1);
+        *f2Out = fabs(f2);
+
+        return 1;
+    }
+
+    double eta (double h, double kx, double ky, double lambda, double T, double x, double y, double t, double phase)
+    {
+        double k = sqrt(kx*kx + ky*ky);
+
+        double b22 = B22(h, k);
+        double b24 = B24(h, k);
+        double b33 = B33(h, k);
+        double b35 = B35(h, k);
+        double b44 = B44(h, k);
+        double b55 = B55(h, k);
+
+        double amp1 = lambda/k;
+        double amp2 = (b22*pow(lambda,2)+b24*pow(lambda,4))/k;
+        double amp3 = (b33*pow(lambda,3)+b35*pow(lambda,5))/k;
+        double amp4 = b44*pow(lambda,4)/k;
+        double amp5 = b55*pow(lambda,5)/k;
+
+        double theta = kx*x + ky*y - 2.0*PII/T*t + phase;
+
+        double C = amp1*cos(theta) 
+		 + amp2*cos(2*theta)
+                 + amp3*cos(3*theta) 
+		 + amp4*cos(4*theta)
+            	 + amp5*cos(5*theta);
+
+        return C;
+    }   
+
+    double U (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z)
+    {
+        double k = sqrt(kx*kx + ky*ky);
+
+        double a11 = A11(d, k);
+        double a13 = A13(d, k);
+        double a15 = A15(d, k);
+        double a22 = A22(d, k);
+        double a24 = A24(d, k);
+        double a33 = A33(d, k);
+        double a35 = A35(d, k);
+        double a44 = A44(d, k);
+        double a55 = A55(d, k);
+
+        double a1u=2.0*PII/T/k*(lambda*a11+pow(lambda,3)*a13+pow(lambda,5)*a15);
+        double a2u=2.0*2.0*PII/T/k*(pow(lambda,2)*a22+pow(lambda,4)*a24);
+        double a3u=3.0*2.0*PII/T/k*(pow(lambda,3)*a33+pow(lambda,5)*a35);
+        double a4u=4.0*2.0*PII/T/k*(pow(lambda,4)*a44);
+        double a5u=5.0*2.0*PII/T/k*(pow(lambda,5)*a55);
+
+        double velU = 0;
+
+        double theta = kx*x + ky*y - 2.0*PII/T*t + phase;
+
+        velU = a1u*cosh(k*z)*cos(theta)
+                 + a2u*cosh(2.0*k*z)*cos(2.0*(theta))
+                 + a3u*cosh(3.0*k*z)*cos(3.0*(theta))
+                 + a4u*cosh(4.0*k*z)*cos(4.0*(theta))
+                 + a5u*cosh(5.0*k*z)*cos(5.0*(theta));
+
+        return velU;
+    }
+
+    double W (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z)
+    {
+        double k = sqrt(kx*kx + ky*ky);
+
+        double a11 = A11(d, k);
+        double a13 = A13(d, k);
+        double a15 = A15(d, k);
+        double a22 = A22(d, k);
+        double a24 = A24(d, k);
+        double a33 = A33(d, k);
+        double a35 = A35(d, k);
+        double a44 = A44(d, k);
+        double a55 = A55(d, k);
+
+        double a1u=2.0*PII/T/k*(lambda*a11+pow(lambda,3)*a13+pow(lambda,5)*a15);
+        double a2u=2.0*2.0*PII/T/k*(pow(lambda,2)*a22+pow(lambda,4)*a24);
+        double a3u=3.0*2.0*PII/T/k*(pow(lambda,3)*a33+pow(lambda,5)*a35);
+        double a4u=4.0*2.0*PII/T/k*(pow(lambda,4)*a44);
+        double a5u=5.0*2.0*PII/T/k*(pow(lambda,5)*a55);
+
+        double velV = 0;
+
+        double theta = kx*x+ky*y-2.0*PII/T*t+phase;
+
+        velV = a1u*sinh(k*z)*sin(theta)
+                + a2u*sinh(2.0*k*z)*sin(2.0*(theta))
+                + a3u*sinh(3.0*k*z)*sin(3.0*(theta))
+                + a4u*sinh(4.0*k*z)*sin(4.0*(theta))
+                + a5u*sinh(5.0*k*z)*sin(5.0*(theta));
+
+        return velV;
+    }
+}
+
+namespace BoussinesqFun
+{
+    #define PII 3.1415926535897932384626433832795028
+    #define G 9.81 
+    
+    double eta (double H, double h, double x, double y, double theta, double t, double X0)
+    {
+        double C = sqrt(G*(H+h));
+        double ts = 3.5*h/sqrt(H/h);
+        double aux = sqrt(3.0*H/(4.0*h))/h;
+        double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta);
+
+        double sup = H * 1.0/pow(cosh( aux * Xa ),2);
+
+        return sup;
+    }
+
+    double Deta1 (double H, double h, double x, double y, double theta, double t, double X0)
+    {
+        double C = sqrt(G*(H+h));
+        double ts = 3.5*h/sqrt(H/h);
+        double aux = sqrt(3.0*H/(4.0*h))/h;
+        double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta);
+
+        double deta = 8.0*aux*h * exp(2.0*aux*Xa) * (1.0-exp(2.0*aux*Xa)) 
+                        / pow(1.0+exp(2.0*aux*Xa),3);
+
+        return deta;
+    }
+
+    double Deta2 (double H, double h, double x, double y, double theta, double t, double X0)
+    {
+        double C = sqrt(G*(H+h));
+        double ts = 3.5*h/sqrt(H/h);
+        double aux = sqrt(3.0*H/(4.0*h))/h;
+        double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta);
+
+        double deta = 16.0*pow(aux,2)*h * exp(2.0*aux*Xa) * (exp(4.0*aux*Xa)
+                        - 4.0*exp(2.0*aux*Xa)+1.0) / pow(1.0+exp(2.0*aux*Xa),4);
+
+        return deta;
+    }
+
+    double Deta3 (double H, double h, double x, double y, double theta, double t, double X0)
+    {
+        double C = sqrt(G*(H+h));
+        double ts = 3.5*h/sqrt(H/h);
+        double aux = sqrt(3.0*H/(4.0*h))/h;
+        double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta);
+
+        double deta = -32.0*pow(aux,3)*h * exp(2.0*aux*Xa) * (exp(6.0*aux*Xa)
+                        - 11.0*exp(4.0*aux*Xa) + 11.0*exp(2.0*aux*Xa)-1.0) / pow(1.0+exp(2.0*aux*Xa),5);
+
+        return deta;
+    }
+
+    double U (double H, double h, double x, double y, double theta, double t, double X0, double z)
+    {
+        double C = sqrt(G*(H+h));
+        double etaSolit = eta( H, h, x, y, theta, t, X0);
+        double Detas2 = Deta2( H, h, x, y, theta, t, X0);
+
+        double vel = C*etaSolit/h
+                *(1.0 - etaSolit/(4.0*h) +
+                pow(h,2)/(3.0*etaSolit)
+                *(1.0 - 3.0/2.0*pow(z/h,2))
+                *Detas2);
+
+        return vel;
+    }
+
+    double W (double H, double h, double x, double y, double theta, double t, double X0, double z)
+    {
+        double C = sqrt(G*(H+h));
+        double etaSolit = eta( H, h, x, y, theta, t, X0);
+        double Detas1 = Deta1( H, h, x, y, theta, t, X0);
+        double Detas3 = Deta3( H, h, x, y, theta, t, X0);
+
+        double vel = -C*z/h
+                *((1.0 - etaSolit/(2.0*h))*Detas1 +
+                pow(h,2)/3.0
+                *(1.0 - 1.0/2.0*pow(z/h,2))
+                *Detas3);
+
+        return vel;
+    }
+}
+

integration/genAbs/common/waveFun.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#ifndef waveFun_H
+#define waveFun_H
+
+namespace otherFun
+{
+    double interpolation (double x1, double x2, double y1, double y2, double xInt);
+}
+
+namespace StokesIFun
+{
+    double waveLength (double h, double T);
+    double eta (double H, double Kx, double x, double Ky, double y, double omega, double t, double phase);
+    double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z);
+    double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z);
+}
+
+namespace StokesIIFun
+{
+    double eta (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase);
+    double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z);
+    double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z);
+    double timeLag (double H, double h, double Kx, double x, double Ky, double y, double T, double phase);
+}
+
+namespace Elliptic
+{
+    void ellipticIntegralsKE (double m, double* K, double* E);
+    void JacobiSnCnDn (double u, double m, double* Sn, double* Cn, double* Dn);
+}
+
+namespace cnoidalFun
+{
+    double eta (double H, double m, double kx, double ky, double T, double x, double y, double t);
+    double etaCnoidal1D (double H, double m, double t, double T, double K, double E);
+    double etaMeanSq (double H, double m, double T);
+    double d1EtaDx (double H, double m, double uCnoidal, double L, double K, double E);
+    double d2EtaDx (double H, double m, double uCnoidal, double L, double K, double E);
+    double d3EtaDx (double H, double m, double uCnoidal, double L, double K, double E);
+    int calculations (double H, double d, double T, double* mOut, double* LOut);
+    double U (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z);
+    double W (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z);
+}
+
+namespace stokesVFun
+{
+    double A11 (double h, double k);
+    double A13 (double h, double k);
+    double A15 (double h, double k);
+    double A22 (double h, double k);
+    double A24 (double h, double k);
+    double A33 (double h, double k);
+    double A35 (double h, double k);
+    double A44 (double h, double k);
+    double A55 (double h, double k);
+    double B22 (double h, double k);
+    double B24 (double h, double k);
+    double B33 (double h, double k);
+    double B33k (double h, double k);
+    double B35 (double h, double k);
+    double B35k (double h, double k);
+    double B44 (double h, double k);
+    double B55 (double h, double k);
+    double B55k (double h, double k);
+    double C1 (double h, double k);
+    double C1k (double h, double k);
+    double C2 (double h, double k);
+    double C2k (double h, double k);
+    double C3 (double h, double k);
+    double C4 (double h, double k);
+    int StokesVNR (double H, double d, double T, double* kOut, double* LambdaOut, double* f1Out, double* f2Out );
+    int StokesExtendedSolver (double H, double d, double T, double* kOut, double* LambdaOut, double* LambdaErrOut );
+    double eta (double h, double kx, double ky, double lambda, double T, double x, double y, double t, double phase);
+    double U (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z);
+    double W (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z);
+}
+
+namespace BoussinesqFun
+{
+    double eta (double H, double h, double x, double y, double theta, double t, double X0);
+    double Deta1 (double H, double h, double x, double y, double theta, double t, double X0);
+    double Deta2 (double H, double h, double x, double y, double theta, double t, double X0);
+    double Deta3 (double H, double h, double x, double y, double theta, double t, double X0);
+    double U (double H, double h, double x, double y, double theta, double t, double X0, double z);
+    double W (double H, double h, double x, double y, double theta, double t, double X0, double z);
+}
+
+#endif

integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C

@@ -0,0 +1,339 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#include "IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(p, iF),
+    nPaddles_(1),
+    leftORright_(-1),
+    waterDepth_(-1),
+    initialDepthABS_(-1),
+    RealwaterDepth_(-1),
+    allCheck_(true),
+    waveDictName_("IHWavesDict")
+{}
+
+
+Foam::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+(
+    const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
+    nPaddles_(ptf.nPaddles_),
+    leftORright_(ptf.leftORright_),    
+    waterDepth_(ptf.waterDepth_),
+    initialDepthABS_(ptf.initialDepthABS_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    allCheck_(ptf.allCheck_),
+    waveDictName_(ptf.waveDictName_)
+{}
+
+
+Foam::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<vector>(p, iF, dict),
+    nPaddles_(dict.lookupOrDefault<label>("nPaddles", 1)),
+    leftORright_(dict.lookupOrDefault<scalar>("leftORright",-1)),
+    waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1 )),
+    initialDepthABS_(dict.lookupOrDefault<scalar>("initialDepthABS",-1)),
+    RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1)),
+    allCheck_(dict.lookupOrDefault<bool>("allCheck", true )),
+    waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict"))
+{}
+
+
+Foam::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+(
+    const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf
+)
+:
+    fixedValueFvPatchField<vector>(ptf),
+    nPaddles_(ptf.nPaddles_),
+    leftORright_(ptf.leftORright_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepthABS_(ptf.initialDepthABS_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    allCheck_(ptf.allCheck_),
+    waveDictName_(ptf.waveDictName_)
+{}
+
+
+Foam::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+(
+    const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(ptf, iF),
+    nPaddles_(ptf.nPaddles_),
+    leftORright_(ptf.leftORright_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepthABS_(ptf.initialDepthABS_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    allCheck_(ptf.allCheck_),
+    waveDictName_(ptf.waveDictName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Variables & constants
+    const vector cMin = gMin(patch().patch().localPoints());
+    const vector cMax = gMax(patch().patch().localPoints());
+    const vector cSpan = cMax - cMin;
+    const scalar zSpan = cSpan[2];
+    scalar dMin = 0.0;
+    scalar dSpan = 0.0;
+    const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan );
+    const volScalarField& alpha = 
+        db().lookupObject<volScalarField>(alphaName());
+    const volVectorField& U = db().lookupObject<volVectorField>("U");
+    const fvMesh& mesh = alpha.mesh();
+    const word& patchName = this->patch().name();
+    const label patchID = mesh.boundaryMesh().findPatchID(patchName);
+    const label nF = patch().faceCells().size();
+    const scalarField alphaCell = 
+        alpha.boundaryField()[patchID].patchInternalField();
+    const vectorField UCell = U.boundaryField()[patchID].patchInternalField();
+    scalarField patchUc = Foam::scalarField(nF, 0.0);
+    scalarField patchVc = Foam::scalarField(nF, 0.0);
+    scalarField patchWc = Foam::scalarField(nF, 0.0);
+    const scalarField patchHeight = patch().Cf().component(2);
+    const scalar g = 9.81;
+
+    // Calculate Z bounds of the faces
+    scalarField zSup, zInf;
+    faceBoundsZ( &zSup, &zInf );
+
+    // Define dictionary
+    IOdictionary IHWavesDict
+    (
+        IOobject
+        (
+            waveDictName_,
+            this->db().time().constant(),
+            this->db(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    scalar currTime = this->db().time().value();
+
+    // Check for errors - Just the first time
+    if (allCheck_)
+    {
+        waterDepth_ = (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1.0));
+
+	initialDepthABS_ = 
+	    (IHWavesDict.lookupOrDefault<scalar>("initialDepthABS",-1.0));
+
+        // Check if the value of nPaddles is correct for the number of columns
+        if (nPaddles_ < 1)
+        {
+            FatalError << "Check nPaddles value." << exit(FatalError);
+        }
+
+        if ( nPaddles_ > 1 )
+        {
+            nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan );
+            reduce(nPaddles_, minOp<label>());
+        }
+    }
+
+    // Grouping part
+    labelList faceGroup = Foam::labelList(nF, 0);
+    scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); 
+    scalarList xGroup = Foam::scalarList(nPaddles_, 0.0);
+    scalarList yGroup = Foam::scalarList(nPaddles_, 0.0);
+
+    for (label i=0; i<nPaddles_; i++)
+    {
+        // Breakpoints, X & Y centre of the paddles
+        dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1);
+        xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) 
+	    + cSpan[0]/(nPaddles_)*i;
+        yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_) 
+	    + cSpan[1]/(nPaddles_)*i;
+    }
+
+    forAll(patchD, patchCells) 
+    {
+        for (label i=0; i<nPaddles_; i++)
+        {
+            if ( (patchD[patchCells]>=dBreakPoints[i])
+                && (patchD[patchCells]<dBreakPoints[i+1]) )
+            {
+                faceGroup[patchCells] = i+1;
+                continue;
+            }
+        }      
+    }
+
+    if (allCheck_)
+    {
+	if (RealwaterDepth_ == -1.0)
+	{
+	    if (waterDepth_ == -1.0)
+	    {
+	        RealwaterDepth_ = 
+		calcWL( alphaCell, faceGroup, zSpan )[0];
+
+		if (initialDepthABS_ !=-1.0)
+		{
+		    RealwaterDepth_ = RealwaterDepth_
+		        + initialDepthABS_;
+		}
+	    }
+	    else if ( waterDepth_ != -1.0 )
+	    {
+	        RealwaterDepth_ = waterDepth_;
+	    }
+	}
+
+        allCheck_ = false;
+    }
+
+    // Calculate water measured levels
+    scalarList measuredLevelsABS = calcWL( alphaCell, faceGroup, zSpan );
+
+    if (initialDepthABS_ !=-1.0)
+    {
+    	forAll(measuredLevelsABS, iterMin)
+    	{
+            measuredLevelsABS[iterMin] = measuredLevelsABS[iterMin] 
+	        + initialDepthABS_;
+    	}
+    }
+
+    forAll(patchHeight, cellIndex)    
+    {
+    	if ( zInf[cellIndex] >= measuredLevelsABS[faceGroup[cellIndex]-1] )
+	{
+	    patchUc[cellIndex] = 0.0;
+	    patchVc[cellIndex] = 0.0;
+             patchWc[cellIndex] = 0.0;
+	}
+	else 
+	{
+	    patchUc[cellIndex] = 
+	        ( -measuredLevelsABS[faceGroup[cellIndex]-1]
+		+ RealwaterDepth_ )
+		* sqrt(g/measuredLevelsABS[faceGroup[cellIndex]-1]);
+
+	    patchUc[cellIndex] = leftORright_ * patchUc[cellIndex];
+            patchVc[cellIndex] = 0.0;
+            patchWc[cellIndex] = 0.0;
+	}
+    }
+
+    const vectorField n1 = Foam::vectorField(nF, vector(1.0, 0.0, 0.0));
+    const vectorField n2 = Foam::vectorField(nF, vector(0.0, 1.0, 0.0));
+    const vectorField n3 = Foam::vectorField(nF, vector(0.0, 0.0, 1.0));
+
+    operator == (n1*patchUc + n2*patchVc + n3*patchWc);
+
+    fixedValueFvPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::
+write(Ostream& os) const
+{
+    fvPatchField<vector>::write(os);
+
+    os.writeKeyword("waveDictName") << waveDictName_ << 
+        token::END_STATEMENT << nl;
+
+    os.writeKeyword("leftORright") << leftORright_ << 
+        token::END_STATEMENT << nl;
+
+    os.writeKeyword("nPaddles") << nPaddles_ << token::END_STATEMENT << nl;
+
+    os.writeKeyword("RealwaterDepth") << RealwaterDepth_ << 
+        token::END_STATEMENT << nl;
+
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //

integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+
+Description
+    Wave absorption boundary condition based on shallow water theory and on a
+    3D approach. Works both in 2D and 3D and for waves out of the shallow water
+    regime.
+
+    Example of the boundary condition specification:
+    @verbatim
+    outlet
+    {
+        type            IH_3D_3DAbsorption_InletVelocityV2;
+        nPaddles        1;  
+        value           uniform (0 0 0);
+	leftORright	-1.0;
+    }
+    @endverbatim
+
+    \heading Function object usage
+    \table
+        Property     | Description             | Required    | Default value
+        type         | type:    IH_3D_3DAbsorption_InletVelocityV2 | yes
+        nPaddles     | number of wavepaddles for absorption | yes | 1  
+	leftORright  | Define location of Boundary condition: Left(1) or Right (-1) | yes | -1
+    \endtable
+
+Note
+    -
+
+SourceFiles
+    IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#ifndef IH_3D_3DAbsorption_InletVelocityFvPatchVectorField_H
+#define IH_3D_3DAbsorption_InletVelocityFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+               Class IH_3D_3DAbsorption_InletVelocityFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Number of paddles
+        label nPaddles_;
+ 
+        //Aborption: Left(1) or Right(-1)
+	scalar leftORright_;
+
+        //- Theoretical Water depth (meters)
+        scalar waterDepth_;        
+
+        //- Numerical Water depth (meters)
+        scalar RealwaterDepth_;
+
+        //- BC check and read values just the first time
+        bool allCheck_;
+
+        //- Initial Depth on BC for absorption (meters)
+        scalar initialDepthABS_;
+
+        //- Dictionary name
+        word waveDictName_;
+
+public:
+
+   //- Runtime type information
+   TypeName("IH_3D_3DAbsorption_InletVelocityV2");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        //  onto a new patch
+        IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        (
+            const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        (
+            const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new IH_3D_3DAbsorption_InletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+        (
+            const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new IH_3D_3DAbsorption_InletVelocityFvPatchVectorField
+                    (*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+    // Other common member functions
+        #include "memberFun.H"
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

integration/genAbs/waveAbsorption/Make/files

@@ -0,0 +1,3 @@
+IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C
+
+LIB = $(FOAM_USER_LIBBIN)/libIHwaveAbsorptionV2esi

integration/genAbs/waveAbsorption/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -DOFVERSION=$(OF_VERSION) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I../common
+
+LIB_LIBS = \
+    -lfiniteVolume

integration/genAbs/waveAbsorption/localMake

@@ -0,0 +1,25 @@
+#!/bin/bash
+
+if [ $WM_PROJECT == "foam" ]; then
+    ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'`
+else
+    ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'`
+fi
+
+#IHC_dbg
+echo $ofversion
+#----
+
+export OF_VERSION=$ofversion
+
+wclean
+
+wmake libso
+
+if (( $? )) ; then 
+    echo "\n\nIH Wave absorption boundary conditions (V2.0) compilation failed"
+    exit; else
+    echo -e "\n\nIH Wave absorption boundary conditions (V2.0) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n";  
+fi
+
+wclean

integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C

@@ -0,0 +1,460 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#include "IH_Waves_InletAlphaFvPatchScalarField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+#include "Random.H"
+
+#include "waveFun.H"
+
+//C++
+#include <stdio.h>
+#include <unistd.h>
+#include <errno.h>
+#include "SquareMatrix.H"
+#include "vector.H"
+#include "Matrix.H"
+
+#include "mpi.h"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::
+IH_Waves_InletAlphaFvPatchScalarField::
+IH_Waves_InletAlphaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(p, iF),
+    wavePeriod_(-1),
+    waveHeight_(-1),
+    waveLength_(-1),
+    waterDepth_(-1),
+    initialDepth_(-1),
+    wavePhase_(3.0*PI()/2.0),
+    lambdaStokesV_(-1),
+    mCnoidal_(-1),
+    nPaddles_(1),
+    tSmooth_(-1),
+    waveDictName_("IHWavesDict"),
+    waveType_("aaa"),
+    waveTheory_("aaa"),
+    allCheck_(true),
+    waveDir_(0),
+    RealwaterDepth_(-1)
+{}
+
+
+Foam::
+IH_Waves_InletAlphaFvPatchScalarField::
+IH_Waves_InletAlphaFvPatchScalarField
+(
+    const IH_Waves_InletAlphaFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_),
+    RealwaterDepth_(ptf.RealwaterDepth_)
+{}
+
+
+Foam::
+IH_Waves_InletAlphaFvPatchScalarField::
+IH_Waves_InletAlphaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<scalar>(p,iF,dict),
+    wavePeriod_(dict.lookupOrDefault<scalar>("wavePeriod",-1)),
+    waveHeight_(dict.lookupOrDefault<scalar>("waveHeight",-1)),
+    waveLength_(dict.lookupOrDefault<scalar>("waveLength",-1)),
+    waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1)),
+    initialDepth_(dict.lookupOrDefault<scalar>("initialDepth",-1)),
+    wavePhase_(dict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)),
+    lambdaStokesV_(dict.lookupOrDefault<scalar>("lambdaStokesV",-1)),
+    mCnoidal_(dict.lookupOrDefault<scalar>("mCnoidal",-1)),
+    nPaddles_(dict.lookupOrDefault<label>("nPaddles",1)),
+    tSmooth_(dict.lookupOrDefault<scalar>("tSmooth",-1)),
+    waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict")),
+    waveType_(dict.lookupOrDefault<word>("waveType","aaa")),
+    waveTheory_(dict.lookupOrDefault<word>("waveTheory","aaa")),
+    allCheck_(dict.lookupOrDefault<bool>("allCheck",true)),
+    waveDir_(dict.lookupOrDefault<scalar>("waveDir",0)),
+    RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1))
+{}
+
+
+Foam::
+IH_Waves_InletAlphaFvPatchScalarField::
+IH_Waves_InletAlphaFvPatchScalarField
+(
+    const IH_Waves_InletAlphaFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchField<scalar>(ptf),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_),
+    RealwaterDepth_(ptf.RealwaterDepth_)
+{}
+
+
+Foam::
+IH_Waves_InletAlphaFvPatchScalarField::
+IH_Waves_InletAlphaFvPatchScalarField
+(
+    const IH_Waves_InletAlphaFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(ptf, iF),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_),
+    RealwaterDepth_(ptf.RealwaterDepth_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::IH_Waves_InletAlphaFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Aux. values
+    scalar auxiliar = 0;
+    scalar auxiliarTotal = 0;
+    scalar X0 = 0;
+    scalarField patchXsolit;
+
+    const vector cMin = gMin(patch().patch().localPoints());
+    const vector cMax = gMax(patch().patch().localPoints());
+    const vector cSpan = cMax - cMin;
+    const scalar zSpan = cSpan[2];
+
+    scalar dMin = 0.0;
+    scalar dSpan = 0.0;
+    const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan );
+
+    // Variables & constants
+    const volScalarField& alpha = 
+        db().lookupObject<volScalarField>(alphaName());
+    const volVectorField& U = db().lookupObject<volVectorField>("U");
+    const fvMesh& mesh = alpha.mesh();
+    const word& patchName = this->patch().name();
+    const label patchID = mesh.boundaryMesh().findPatchID(patchName);
+    const label nF = patch().faceCells().size();
+    labelList cellGroup = Foam::labelList(nF, 1);
+    const scalarField alphaCell = 
+        alpha.boundaryField()[patchID].patchInternalField();
+    const vectorField UCell = U.boundaryField()[patchID].patchInternalField();
+    scalarField patchVOF = alphaCell; 
+    const labelList celdas = patch().faceCells();
+    const scalarField patchHeight = patch().Cf().component(2);
+
+    // Calculate Z bounds of the faces
+    scalarField zSup, zInf;
+    faceBoundsZ( &zSup, &zInf );
+
+    // Wave variables
+    scalar waveOmega;
+    scalarList waveOmegas;
+    scalar waveK;
+    scalarList waveKs;
+    scalar waveAngle;
+    scalar waveKx;
+    scalarList waveKxs;
+    scalar waveKy;
+
+    // Define dictionary
+    IOdictionary IHWavesDict
+    (
+        IOobject
+        (
+            waveDictName_,
+            this->db().time().constant(),
+            this->db(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    scalar currTime = this->db().time().value();
+
+    // Check for errors - Just the first time
+    if (allCheck_)
+    {
+        waveType_ = (IHWavesDict.lookupOrDefault<word>("waveType","aaa")); 
+        tSmooth_ = (IHWavesDict.lookupOrDefault<scalar>("tSmooth",-1.0));
+	waterDepth_ = 
+	    (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1.0));
+	initialDepth_ = 
+            (IHWavesDict.lookupOrDefault<scalar>("initialDepth",-1.0));
+    	nPaddles_ = (IHWavesDict.lookupOrDefault<label>("nPaddles",1));
+    } 
+
+    // Grouping part
+    scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); 
+    scalarList xGroup = Foam::scalarList(nPaddles_, 0.0);
+    scalarList yGroup = Foam::scalarList(nPaddles_, 0.0);
+
+    for (label i=0; i<nPaddles_; i++)
+    {
+        // Breakpoints, X & Y centre of the paddles
+        dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1);
+        xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) 
+	    + cSpan[0]/(nPaddles_)*i;
+        yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_)
+	    + cSpan[1]/(nPaddles_)*i;
+    }
+
+    forAll(patchD, patchCells) 
+    {
+        for (label i=0; i<nPaddles_; i++)
+        {
+            if ( (patchD[patchCells]>=dBreakPoints[i])
+	        && (patchD[patchCells]<dBreakPoints[i+1]) )
+            {
+                cellGroup[patchCells] = i+1;
+                continue;
+            }
+        }      
+    }
+
+    // Check for errors - Just the first time
+    if (allCheck_)
+    {
+	if (RealwaterDepth_ == -1.0)
+	{
+	    if (waterDepth_ == -1.0)
+	    {
+	        RealwaterDepth_ = 
+		    calcWL( alphaCell, cellGroup, zSpan )[0];
+		if (initialDepth_ !=-1.0)
+		{
+			RealwaterDepth_ = RealwaterDepth_
+			    + initialDepth_;
+		}
+	    }
+	    else if ( waterDepth_ != -1.0 )
+	    {
+		RealwaterDepth_ = waterDepth_;
+	    }
+	}
+    }
+
+    if (allCheck_)
+    {
+	waveTheory_ = (IHWavesDict.lookupOrDefault<word>("waveTheory","aaa"));
+    	waveHeight_ = (IHWavesDict.lookupOrDefault<scalar>("waveHeight",-1));
+    	wavePeriod_ = (IHWavesDict.lookupOrDefault<scalar>("wavePeriod",-1));
+    	waveDir_ = (IHWavesDict.lookupOrDefault<scalar>("waveDir",0));
+    	wavePhase_ = 
+	    (IHWavesDict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0));
+	
+	if ( waveType_ == "regular" )
+	{
+	   waveLength_ = StokesIFun::waveLength (RealwaterDepth_, wavePeriod_);
+	}
+    }
+
+    if ( waveType_ == "regular" )
+    {
+        waveOmega = (2.0*PI())/wavePeriod_;
+        waveK = 2.0*PI()/waveLength_;
+
+        waveAngle = waveDir_*PI()/180.0;
+        waveKx = waveK*cos(waveAngle);
+        waveKy = waveK*sin(waveAngle);
+    }
+    else if ( waveType_ == "solitary" )
+    {
+        waveAngle = waveDir_*PI()/180.0;
+        patchXsolit = patch().Cf().component(0)*cos(waveAngle) 
+	    + patch().Cf().component(1)*sin(waveAngle);
+        X0 = gMin(patchXsolit);
+    }
+
+    scalar timeMult = 1.0;
+
+    if ( tSmooth_ > 0 && currTime < tSmooth_ )
+    {
+	timeMult = currTime/tSmooth_;
+    }
+
+    if (allCheck_)
+    {
+        if ( waveType_ == "regular" )
+        {
+            #include "checkInputErrorsRegular.H"
+        }
+        else if ( waveType_ == "solitary" )
+        {
+            #include "checkInputErrorsSolitary.H"
+        }
+        else
+        {
+            FatalError << "Wave type not supported, use:\n" 
+	        << "regular, solitary" << exit(FatalError);
+        }
+	allCheck_ = false; 
+    }
+
+    scalarList calculatedLevel (nPaddles_,0.0);
+
+    if ( waveType_ == "regular" )
+    {
+	if (waveDir_ == 0)
+	{
+           #include "calculatedLevelRegularNormal.H"
+	}
+	else
+	{
+	   #include "calculatedLevelRegular.H"
+	}
+    }
+    else if ( waveType_ == "solitary" )
+    {
+        #include "calculatedLevelSolitary.H"
+    }
+
+    forAll(patchHeight, cellIndex)    
+    {
+	auxiliarTotal = 0;
+	auxiliar = 0;
+	
+	if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+	{
+	    patchVOF[cellIndex] = 1.0;
+	}
+	else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+	{
+	    patchVOF[cellIndex] = 0.0;
+	}
+	else 
+	{
+	    auxiliar = calculatedLevel[cellGroup[cellIndex]-1]
+	        - zInf[cellIndex];
+	    auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+	    auxiliarTotal = auxiliar/auxiliarTotal;
+	    patchVOF[cellIndex] = auxiliarTotal;
+	}
+    }
+
+    operator == (patchVOF);
+
+    fixedValueFvPatchField<scalar>::updateCoeffs();
+}
+
+
+void Foam::IH_Waves_InletAlphaFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+
+    os.writeKeyword("waveDictName") << waveDictName_ << 
+        token::END_STATEMENT << nl;
+
+    os.writeKeyword("RealwaterDepth") << RealwaterDepth_ << 
+        token::END_STATEMENT << nl;
+
+    os.writeKeyword("initialDepth") << initialDepth_ << 
+        token::END_STATEMENT << nl;
+
+    writeEntry("value", os);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchScalarField,
+       IH_Waves_InletAlphaFvPatchScalarField
+   );
+
+}
+
+
+// ************************************************************************* //

integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.H

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::IH_Waves_InletAlphaFvPatchScalarField
+
+Description
+    Describes a volumetric/mass flow normal Scalar boundary condition by its
+    magnitude as an integral over its area.
+
+    The basis of the patch (volumetric or mass) is determined by the
+    dimensions of the flux, phi.
+    The current density is used to correct the Alpha when applying the
+    mass basis.
+
+    Example of the boundary condition specification:
+    @verbatim
+    inlet
+    {
+        type            IH_Waves_InletAlphaV2;
+        waveDict        IHWavesDict;
+    }
+    @endverbatim
+
+    \heading Function object usage
+    \table
+        Property     | Description             | Required    | Default value
+        type         | type:    IH_Waves_InletAlphaV2 | yes
+        waveDict     | Dictionary where variables for generation/absorption are defined | yes | IHWavesDict  
+    \endtable
+
+Note
+    - 
+
+SourceFiles
+    IH_Waves_InletAlphaFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#ifndef IH_Waves_InletAlphaFvPatchScalarField_H
+#define IH_Waves_InletAlphaFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+               Class IH_Waves_InletAlphaFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class IH_Waves_InletAlphaFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Wave period (seconds)
+        scalar wavePeriod_;
+
+        //- Wave height (meters)
+        scalar waveHeight_;
+
+        //- Wave length (meters)
+        scalar waveLength_;
+
+        //- Water depth (meters)
+        scalar waterDepth_;
+
+        //- Wave phase (radians)
+        scalar wavePhase_;
+
+        //- Lambda StokesV parameter
+        scalar lambdaStokesV_;
+
+        //- m Cnoidal parameter
+        scalar mCnoidal_;
+
+        //- Number of different paddles (for absorption)
+        label nPaddles_;
+
+        //- Smoothing factor for the wave(s) (seconds)
+        scalar tSmooth_;
+
+        //- Dictionary name
+        word waveDictName_;
+
+        //- Regular or Irregular wave(s)
+        word waveType_;
+
+        //- Name of the theory
+        word waveTheory_;
+
+        //- BC check and read values just the first time
+        bool allCheck_;
+
+        //- Direction of propagation (degrees, from X axis)
+        scalar waveDir_;
+
+        //- Numerical Water depth (meters)
+	scalar RealwaterDepth_;
+
+        //- Initial Depth on BC for wave generation (meters)
+        scalar initialDepth_;
+
+public:
+
+   //- Runtime type information
+   TypeName("IH_Waves_InletAlphaV2");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        IH_Waves_InletAlphaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        IH_Waves_InletAlphaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  IH_Waves_InletAlphaFvPatchScalarField
+        //  onto a new patch
+        IH_Waves_InletAlphaFvPatchScalarField
+        (
+            const IH_Waves_InletAlphaFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        IH_Waves_InletAlphaFvPatchScalarField
+        (
+            const IH_Waves_InletAlphaFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new IH_Waves_InletAlphaFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        IH_Waves_InletAlphaFvPatchScalarField
+        (
+            const IH_Waves_InletAlphaFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new IH_Waves_InletAlphaFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+    // Other common member functions
+        #include "memberFun.H"
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C

@@ -0,0 +1,495 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#include "IH_Waves_InletVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+
+#include "waveFun.H"
+
+#include <stdio.h>
+#include <unistd.h>
+#include <errno.h>
+#include "SquareMatrix.H"
+#include "vector.H"
+#include "Matrix.H"
+
+#include "mpi.h"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::
+IH_Waves_InletVelocityFvPatchVectorField::
+IH_Waves_InletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(p, iF),
+    wavePeriod_(-1),
+    waveHeight_(-1),
+    waveLength_(-1),
+    waterDepth_(-1),
+    initialDepth_(-1),
+    RealwaterDepth_(-1),
+    wavePhase_(3.0*PI()/2.0),
+    lambdaStokesV_(-1),
+    mCnoidal_(-1),
+    genAbs_(false),
+    nPaddles_(1),
+    tSmooth_(-1),
+    leftORright_(-1),
+    waveDictName_("IHWavesDict"),
+    waveType_("aaa"),
+    waveTheory_("aaa"),
+    allCheck_(true),
+    waveDir_(0)
+{}
+
+
+Foam::
+IH_Waves_InletVelocityFvPatchVectorField::
+IH_Waves_InletVelocityFvPatchVectorField
+(
+    const IH_Waves_InletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    genAbs_(ptf.genAbs_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    leftORright_(ptf.leftORright_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_)
+{}
+
+
+Foam::
+IH_Waves_InletVelocityFvPatchVectorField::
+IH_Waves_InletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<vector>(p, iF, dict),
+    wavePeriod_(dict.lookupOrDefault<scalar>("wavePeriod",-1)),
+    waveHeight_(dict.lookupOrDefault<scalar>("waveHeight",-1)),
+    waveLength_(dict.lookupOrDefault<scalar>("waveLength",-1)),
+    waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1)),
+    initialDepth_(dict.lookupOrDefault<scalar>("initialDepth",-1)),
+    RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1)),
+    wavePhase_(dict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)),
+    lambdaStokesV_(dict.lookupOrDefault<scalar>("lambdaStokesV",-1)),
+    mCnoidal_(dict.lookupOrDefault<scalar>("mCnoidal",-1)),
+    genAbs_(dict.lookupOrDefault<bool>("genAbs",false)),
+    nPaddles_(dict.lookupOrDefault<label>("nPaddles",1)),
+    tSmooth_(dict.lookupOrDefault<scalar>("tSmooth",-1)),
+    leftORright_(dict.lookupOrDefault<scalar>("leftORright",-1)),
+    waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict")),
+    waveType_(dict.lookupOrDefault<word>("waveType","aaa")),
+    waveTheory_(dict.lookupOrDefault<word>("waveTheory","aaa")),
+    allCheck_(dict.lookupOrDefault<bool>("allCheck",true)),
+    waveDir_(dict.lookupOrDefault<scalar>("waveDir",0))
+{}
+
+
+Foam::
+IH_Waves_InletVelocityFvPatchVectorField::
+IH_Waves_InletVelocityFvPatchVectorField
+(
+    const IH_Waves_InletVelocityFvPatchVectorField& ptf
+)
+:
+    fixedValueFvPatchField<vector>(ptf),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    genAbs_(ptf.genAbs_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    leftORright_(ptf.leftORright_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_)
+{}
+
+
+Foam::
+IH_Waves_InletVelocityFvPatchVectorField::
+IH_Waves_InletVelocityFvPatchVectorField
+(
+    const IH_Waves_InletVelocityFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(ptf, iF),
+    wavePeriod_(ptf.wavePeriod_),
+    waveHeight_(ptf.waveHeight_),
+    waveLength_(ptf.waveLength_),
+    waterDepth_(ptf.waterDepth_),
+    initialDepth_(ptf.initialDepth_),
+    RealwaterDepth_(ptf.RealwaterDepth_),
+    wavePhase_(ptf.wavePhase_),
+    lambdaStokesV_(ptf.lambdaStokesV_),
+    mCnoidal_(ptf.mCnoidal_),
+    genAbs_(ptf.genAbs_),
+    nPaddles_(ptf.nPaddles_),
+    tSmooth_(ptf.tSmooth_),
+    leftORright_(ptf.leftORright_),
+    waveDictName_(ptf.waveDictName_),
+    waveType_(ptf.waveType_),
+    waveTheory_(ptf.waveTheory_),
+    allCheck_(ptf.allCheck_),
+    waveDir_(ptf.waveDir_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::IH_Waves_InletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    } 
+
+    // Auxiliar variables
+    scalar auxiliar = 0; 
+    scalar auxiliarTotal = 0;
+    scalar X0 = 0;
+    scalarField patchXsolit;
+
+    // Variables stream function
+    scalar celerity = 0;
+
+    // 3D Variables
+    const vector cMin = gMin(patch().patch().localPoints());
+    const vector cMax = gMax(patch().patch().localPoints());
+    const vector cSpan = cMax - cMin;
+    const scalar zSpan = cSpan[2];
+
+    scalar dMin = 0.0;
+    scalar dSpan = 0.0;
+    const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan );
+
+    // Variables & constants
+    const volScalarField& alpha = 
+        db().lookupObject<volScalarField>(alphaName());
+    const fvMesh& mesh = alpha.mesh();
+    const word& patchName = this->patch().name();
+    const label patchID = mesh.boundaryMesh().findPatchID(patchName);
+    const label nF = patch().faceCells().size();
+    labelList cellGroup = Foam::labelList(nF, 1);
+    const scalarField alphaCell = 
+        alpha.boundaryField()[patchID].patchInternalField();
+    scalarField patchU = Foam::scalarField(nF, 0.0);
+    scalarField patchUABS = Foam::scalarField(nF, 0.0);
+    scalarField patchV = Foam::scalarField(nF, 0.0);
+    scalarField patchVABS = Foam::scalarField(nF, 0.0);
+    scalarField patchW = Foam::scalarField(nF, 0.0);
+    const labelList celdas = patch().faceCells();
+    const scalarField patchHeight = patch().Cf().component(2);
+    const scalar g = 9.81;
+
+    // Calculate Z bounds of the faces
+    scalarField zSup, zInf;
+    faceBoundsZ( &zSup, &zInf );
+
+    // Waves variables
+    scalar waveOmega;
+    scalar waveK;
+    scalar waveAngle;
+    scalar waveKx;
+    scalar waveKy;
+
+    // Define dictionary
+    IOdictionary IHWavesDict
+    (
+        IOobject
+        (
+            waveDictName_,
+            this->db().time().constant(),
+            this->db(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    scalar currTime = this->db().time().value();
+
+    // Check for errors - Just the first time
+    if (allCheck_)
+    {
+        waveType_ = (IHWavesDict.lookupOrDefault<word>("waveType","aaa")); 
+        tSmooth_ = (IHWavesDict.lookupOrDefault<scalar>("tSmooth",-1.0));
+	genAbs_ = (IHWavesDict.lookupOrDefault<bool>("genAbs",false));
+	waterDepth_ = (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1));
+	initialDepth_ = 
+	    (IHWavesDict.lookupOrDefault<scalar>("initialDepth",-1));
+    	nPaddles_ = (IHWavesDict.lookupOrDefault<label>("nPaddles",1));
+    }
+
+    // Grouping part
+    scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); 
+    scalarList xGroup = Foam::scalarList(nPaddles_, 0.0);
+    scalarList yGroup = Foam::scalarList(nPaddles_, 0.0);
+
+    for (label i=0; i<nPaddles_; i++)
+    {
+        // Breakpoints, X & Y centre of the paddles
+        dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1);
+        xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) 
+	    + cSpan[0]/(nPaddles_)*i;
+        yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_) 
+	    + cSpan[1]/(nPaddles_)*i;
+    }
+
+    forAll(patchD, patchCells) 
+    {
+        for (label i=0; i<nPaddles_; i++)
+        {
+            if ( (patchD[patchCells]>=dBreakPoints[i]) 
+                && (patchD[patchCells]<dBreakPoints[i+1]) )
+            {
+                cellGroup[patchCells] = i+1; // Group of each face
+                continue;
+            }
+        }      
+    }
+
+    // Check for errors - Just the first time
+    if (allCheck_)
+    {
+	if (RealwaterDepth_ == -1.0)
+	{
+	    if (waterDepth_ == -1.0)
+	    {
+	        RealwaterDepth_ = 
+		    calcWL( alphaCell, cellGroup, zSpan )[0];
+		if (initialDepth_ !=-1.0)
+		    {
+		        RealwaterDepth_ = RealwaterDepth_ 
+			    + initialDepth_;
+		    }
+	    }
+	    else if ( waterDepth_ != -1.0 )
+	    {
+		RealwaterDepth_ = waterDepth_;	
+	    }
+	}
+    }
+
+    if (allCheck_)
+    {
+	waveTheory_ = (IHWavesDict.lookupOrDefault<word>("waveTheory","aaa"));
+    	waveHeight_ = (IHWavesDict.lookupOrDefault<scalar>("waveHeight",-1));
+    	wavePeriod_ = (IHWavesDict.lookupOrDefault<scalar>("wavePeriod",-1));
+    	waveDir_ = (IHWavesDict.lookupOrDefault<scalar>("waveDir",0));
+    	genAbs_ = (IHWavesDict.lookupOrDefault<bool>("genAbs",false));
+    	wavePhase_ = 
+	    (IHWavesDict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0));
+	
+	if ( waveType_ == "regular" )
+	{
+		waveLength_ = StokesIFun::waveLength (
+		    RealwaterDepth_, wavePeriod_ );
+	}
+    }
+
+    if ( waveType_ == "regular" )
+    {
+        waveOmega = (2.0*PI())/wavePeriod_;
+        waveK = 2.0*PI()/waveLength_;
+
+        celerity = waveLength_/wavePeriod_;
+
+        waveAngle = waveDir_*PI()/180.0;
+        waveKx = waveK*cos(waveAngle);
+        waveKy = waveK*sin(waveAngle);
+    }
+    else if ( waveType_ == "solitary" )
+    {
+        waveAngle = waveDir_*PI()/180.0;
+        patchXsolit = patch().Cf().component(0)*cos(waveAngle) 
+            	      + patch().Cf().component(1)*sin(waveAngle);
+        X0 = gMin(patchXsolit);
+    }
+
+    scalar timeMult = 1.0;
+
+    if ( tSmooth_ > 0 && currTime < tSmooth_ )
+    {
+	timeMult = currTime/tSmooth_;
+    }
+
+    if (allCheck_)
+    {
+        if ( waveType_ == "regular" )
+        {
+            #include "checkInputErrorsRegular.H"
+        }
+        else if ( waveType_ == "solitary" )
+        {
+            #include "checkInputErrorsSolitary.H"
+        }
+        else
+        {
+            FatalError << "Wave type not supported, use:\n" 
+	        << "regular, solitary."<< exit(FatalError);
+        }
+
+	allCheck_ = false; 
+    }
+
+    scalarList calculatedLevel (nPaddles_,0.0);  
+
+    if ( waveType_ == "regular" )
+    {
+	if (waveDir_ == 0)
+	{
+           #include "calculatedLevelRegularNormal.H"
+	}
+	else
+	{
+	   #include "calculatedLevelRegular.H"
+	}
+    }
+    else if ( waveType_ == "solitary" )
+    {
+        #include "calculatedLevelSolitary.H"
+    }
+
+    scalarList measuredLevels (nPaddles_,0.0);
+    forAll(measuredLevels, iterMin)
+    {
+        measuredLevels[iterMin] = RealwaterDepth_;
+    }
+
+    scalarList measuredLevelsGENAB = calcWL( alphaCell, cellGroup, zSpan );
+
+    if (initialDepth_ !=-1.0)
+    {
+    	forAll(measuredLevelsGENAB, iterMin)
+    	{
+            measuredLevelsGENAB[iterMin] = measuredLevelsGENAB[iterMin] 
+	        + initialDepth_;
+    	}
+    }
+
+    // Define heights as minimum of calculatedLevel and measuredLevels
+    scalarList heights (nPaddles_,0.0);
+    forAll(heights, iterMin)
+    {
+        heights[iterMin] = 
+	    min(calculatedLevel[iterMin],measuredLevels[iterMin]);
+    }
+    
+    forAll(patchHeight, cellIndex)    
+    {
+    	#include "velocityProfile.H"
+    
+    	patchU[cellIndex] = patchU[cellIndex]*alphaCell[cellIndex];
+    	patchV[cellIndex] = patchV[cellIndex]*alphaCell[cellIndex];
+    	patchW[cellIndex] = patchW[cellIndex]*alphaCell[cellIndex];
+    }
+
+    const vectorField n1 = Foam::vectorField(nF, vector(1.0, 0.0, 0.0));
+    const vectorField n2 = Foam::vectorField(nF, vector(0.0, 1.0, 0.0));
+    const vectorField n3 = Foam::vectorField(nF, vector(0.0, 0.0, 1.0));
+
+    operator == (n1*patchU + n2*patchV + n3*patchW);
+
+    fixedValueFvPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::IH_Waves_InletVelocityFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchField<vector>::write(os);
+
+    os.writeKeyword("waveDictName") << waveDictName_ 
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("leftORright") << leftORright_ 
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("RealwaterDepth") << RealwaterDepth_ 
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("initialDepth") << initialDepth_ 
+        << token::END_STATEMENT << nl;
+
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+   makePatchTypeField
+   (
+       fvPatchVectorField,
+       IH_Waves_InletVelocityFvPatchVectorField
+   );
+}
+
+
+// ************************************************************************* //

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.H

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::IH_Waves_InletVelocityFvPatchVectorField
+
+Description
+    Describes a volumetric/mass flow normal vector boundary condition by its
+    magnitude as an integral over its area.
+
+    The basis of the patch (volumetric or mass) is determined by the
+    dimensions of the flux, phi.
+    The current density is used to correct the velocity when applying the
+    mass basis.
+
+    Example of the boundary condition specification:
+    @verbatim
+    inlet
+    {
+        type            IH_Waves_InletVelocityV2;
+        waveDict        IHWavesDict;
+        value           uniform (0 0 0);
+	leftORright	1.0;
+
+    }
+    @endverbatim
+
+    \heading Function object usage
+    \table
+        Property     | Description             | Required    | Default value
+        type         | type:    IH_Waves_InletVelocityV2 | yes
+        waveDict     | Dictionary where variables for generation/absorption are defined | yes | IHWavesDict 
+	leftORright  | Define location of Boundary condition: Left(1) or Right (-1) | yes | -1 
+    \endtable
+
+Note
+    - The value is positive inwards
+    - May not work correctly for transonic inlets
+    - Strange behaviour with potentialFoam since the U equation is not solved
+
+SourceFiles
+    IH_Waves_InletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+#ifndef IH_Waves_InletVelocityFvPatchVectorField_H
+#define IH_Waves_InletVelocityFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+               Class IH_Waves_InletVelocityFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class IH_Waves_InletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Wave period (seconds)
+        scalar wavePeriod_;
+
+        //- Wave height (meters)
+        scalar waveHeight_;
+
+        //- Wave length (meters)
+        scalar waveLength_;
+
+        //- Water depth (meters)
+        scalar waterDepth_;
+
+        //- Wave phase (radians)
+        scalar wavePhase_;
+
+        //- Lambda StokesV parameter
+        scalar lambdaStokesV_;
+
+        //- m Cnoidal parameter
+        scalar mCnoidal_;
+
+        //- Generation + Absorption at the same time in the boundary (on=1, off=0)
+        bool genAbs_;
+
+        //- Number of different paddles (for absorption)
+        label nPaddles_;
+
+        //- Smoothing factor for the wave(s) (seconds)
+        scalar tSmooth_;
+
+        //Aborption Left(1) or Right(-1)
+	scalar leftORright_;
+
+        //- Dictionary name
+        word waveDictName_;
+
+        //- Regular or Irregular
+        word waveType_;
+
+        //- Name of the theory
+        word waveTheory_;
+
+        //- BC check and read values just the first time
+        bool allCheck_;
+
+        //- Direction of propagation (degrees, from X axis)
+        scalar waveDir_;
+
+        //- Initial water depth (referece)
+	scalar RealwaterDepth_;
+
+        //- Water depth (meters)
+        scalar initialDepth_;
+
+public:
+
+   //- Runtime type information
+   TypeName("IH_Waves_InletVelocityV2");
+
+
+   // Constructors
+
+        //- Construct from patch and internal field
+        IH_Waves_InletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        IH_Waves_InletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  IH_Waves_InletVelocityFvPatchVectorField
+        //  onto a new patch
+        IH_Waves_InletVelocityFvPatchVectorField
+        (
+            const IH_Waves_InletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        IH_Waves_InletVelocityFvPatchVectorField
+        (
+            const IH_Waves_InletVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new IH_Waves_InletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        IH_Waves_InletVelocityFvPatchVectorField
+        (
+            const IH_Waves_InletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new IH_Waves_InletVelocityFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+    // Other common member functions
+        #include "memberFun.H"
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileBoussinesq.H

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = BoussinesqFun :: U 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                xGroup[cellGroup[cellIndex]-1], 
+                yGroup[cellGroup[cellIndex]-1], 
+                waveAngle, 
+                currTime, 
+                X0, 
+                patchHeight[cellIndex]
+             );
+
+        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult;
+        patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult;
+
+        patchW[cellIndex] = BoussinesqFun :: W 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                xGroup[cellGroup[cellIndex]-1], 
+                yGroup[cellGroup[cellIndex]-1], 
+                waveAngle, 
+                currTime, 
+                X0, 
+                patchHeight[cellIndex]
+            );
+
+	patchW[cellIndex] = patchW[cellIndex] * timeMult;
+
+    }
+    else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = 0.0;
+        patchV[cellIndex] = 0.0;
+        patchW[cellIndex] = 0.0;
+    } 
+    else 
+    {	
+        if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    < measuredLevels[cellGroup[cellIndex]-1] )
+        {
+	    if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]))
+	    {
+                auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		    - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = BoussinesqFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveAngle, 
+                    	currTime, 
+                    	X0, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle) * timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle) * timeMult;
+
+		patchW[cellIndex] = BoussinesqFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveAngle, 
+                    	currTime, 
+                    	X0, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex]
+		    * timeMult;	
+	    }
+	}	
+	else if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    > measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+	    {
+	        patchU[cellIndex] = BoussinesqFun :: U 
+		    (
+                        waveHeight_, 
+                        RealwaterDepth_, 
+                        xGroup[cellGroup[cellIndex]-1], 
+                        yGroup[cellGroup[cellIndex]-1], 
+                        waveAngle, 
+                        currTime, 
+                        X0, 
+                        measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) 
+		    * timeMult;
+		patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) 
+		    * timeMult;
+
+		patchW[cellIndex] = BoussinesqFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveAngle, 
+                    	currTime, 
+                    	X0, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchW[cellIndex] = patchW[cellIndex]*timeMult;			
+	    }
+	    else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		     - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = BoussinesqFun :: U 
+		    (
+		        waveHeight_, 
+			RealwaterDepth_, 
+			xGroup[cellGroup[cellIndex]-1], 
+			yGroup[cellGroup[cellIndex]-1], 
+			waveAngle, 
+			currTime, 
+			X0, 
+			measuredLevels[cellGroup[cellIndex]-1]
+		    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle) * timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle) * timeMult;
+
+		patchW[cellIndex] = BoussinesqFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveAngle, 
+                    	currTime, 
+                    	X0, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                   );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+            }
+        }	
+    }

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileCnoidal.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = cnoidalFun :: U 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                mCnoidal_, 
+                waveKx, 
+                waveKy, 
+                wavePeriod_, 
+                xGroup[cellGroup[cellIndex]-1], 
+                yGroup[cellGroup[cellIndex]-1], 
+                currTime, 
+                patchHeight[cellIndex]
+            );
+
+        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult;
+        patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult;
+
+        patchW[cellIndex] = cnoidalFun :: W 
+	    (
+	        waveHeight_, 
+                RealwaterDepth_, 
+                mCnoidal_, 
+                waveKx, 
+                waveKy, 
+                wavePeriod_, 
+                xGroup[cellGroup[cellIndex]-1], 
+                yGroup[cellGroup[cellIndex]-1], 
+                currTime, 
+                patchHeight[cellIndex]
+            );
+
+	patchW[cellIndex] = patchW[cellIndex]*timeMult;
+
+    }
+    else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = 0.0;
+        patchV[cellIndex] = 0.0;
+        patchW[cellIndex] = 0.0;
+    } 
+    else 
+    {	
+        if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    < measuredLevels[cellGroup[cellIndex]-1] )
+        {
+	    if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]))
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = cnoidalFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+                    	waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle)*timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle)*timeMult;
+
+		patchW[cellIndex] = cnoidalFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+			waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+	    }
+	}	
+	else if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    > measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+	    {
+	        patchU[cellIndex] = cnoidalFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+                    	waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) 
+		    * timeMult;
+		patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) 
+		    * timeMult;
+
+		patchW[cellIndex] = cnoidalFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+			waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchW[cellIndex] = patchW[cellIndex]*timeMult;			
+	    }
+	    else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = cnoidalFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+                    	waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle) * timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle) * timeMult;
+
+		patchW[cellIndex] = cnoidalFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	mCnoidal_, 
+			waveKx, 
+                    	waveKy, 
+                    	wavePeriod_, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	currTime, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+	    }
+        }	
+    }

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileGENAB.H

@@ -0,0 +1,20 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if ( genAbs_ )
+    {  
+	if ( zInf[cellIndex] >= measuredLevelsGENAB[cellGroup[cellIndex]-1] )
+	{
+	    patchUABS[cellIndex] = 0.0;
+	}
+	else 
+	{
+	    patchUABS[cellIndex] = ( calculatedLevel[cellGroup[cellIndex]-1] - measuredLevelsGENAB[cellGroup[cellIndex]-1] ) * sqrt(g/measuredLevelsGENAB[cellGroup[cellIndex]-1]);
+        }
+    }
+    patchU[cellIndex] = patchU[cellIndex] + leftORright_ * patchUABS[cellIndex];

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesI.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = StokesIFun :: U 
+	    (
+                 waveHeight_, 
+                 RealwaterDepth_, 
+                 waveKx, 
+                 xGroup[cellGroup[cellIndex]-1], 
+                 waveKy, 
+                 yGroup[cellGroup[cellIndex]-1], 
+                 waveOmega, 
+                 currTime, 
+                 wavePhase_, 
+                 patchHeight[cellIndex]
+            );
+
+        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult;
+        patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult;
+
+        patchW[cellIndex] = StokesIFun :: W 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                waveKx, 
+                xGroup[cellGroup[cellIndex]-1], 
+                waveKy, 
+                yGroup[cellGroup[cellIndex]-1], 
+                waveOmega, 
+                currTime, 
+                wavePhase_, 
+                patchHeight[cellIndex]
+            );
+
+	patchW[cellIndex] = patchW[cellIndex]*timeMult;
+    }
+    else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = 0.0;
+        patchV[cellIndex] = 0.0;
+        patchW[cellIndex] = 0.0;
+    } 
+    else 
+    {	
+        if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    < measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		     - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = StokesIFun :: U 
+		    (
+		        waveHeight_, 
+		        RealwaterDepth_, 
+		        waveKx, 
+		        xGroup[cellGroup[cellIndex]-1], 
+		        waveKy, 
+		        yGroup[cellGroup[cellIndex]-1], 
+		        waveOmega, 
+		        currTime, 
+		        wavePhase_, 
+		        patchHeight[cellIndex]
+		    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle)*timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle)*timeMult;
+
+		patchW[cellIndex] = StokesIFun :: W 
+		    (
+		        waveHeight_, 
+		        RealwaterDepth_, 
+		        waveKx, 
+		        xGroup[cellGroup[cellIndex]-1], 
+		        waveKy, 
+		        yGroup[cellGroup[cellIndex]-1], 
+		        waveOmega, 
+		        currTime, 
+		        wavePhase_, 
+		        patchHeight[cellIndex]
+		    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex]
+		    * timeMult;	
+		
+	    }
+	}
+	else if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    > measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+	    {
+	        patchU[cellIndex] = StokesIFun :: U 
+		    (
+		        waveHeight_, 
+			RealwaterDepth_, 
+			waveKx, 
+			xGroup[cellGroup[cellIndex]-1], 
+			waveKy, 
+			yGroup[cellGroup[cellIndex]-1], 
+			waveOmega, 
+			currTime, 
+			wavePhase_, 
+			measuredLevels[cellGroup[cellIndex]-1]
+		    );
+
+		patchV[cellIndex] = patchU[cellIndex] 
+		    * sin(waveAngle)*timeMult;
+		patchU[cellIndex] = patchU[cellIndex] 
+		    * cos(waveAngle)*timeMult;
+
+		patchW[cellIndex] = StokesIFun :: W 
+		    (
+		        waveHeight_, 
+			RealwaterDepth_, 
+			waveKx, 
+			xGroup[cellGroup[cellIndex]-1], 
+			waveKy, 
+			yGroup[cellGroup[cellIndex]-1], 
+			waveOmega, 
+			currTime, 
+			wavePhase_, 
+			measuredLevels[cellGroup[cellIndex]-1]
+		    );
+
+		patchW[cellIndex] = patchW[cellIndex]*timeMult;	
+		
+	    }
+	    else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]))
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		    - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = StokesIFun :: U 
+		    (
+		        waveHeight_, 
+		        RealwaterDepth_, 
+		        waveKx, 
+		        xGroup[cellGroup[cellIndex]-1], 
+		        waveKy, 
+		        yGroup[cellGroup[cellIndex]-1], 
+		        waveOmega, 
+		        currTime, 
+		        wavePhase_, 
+		        measuredLevels[cellGroup[cellIndex]-1]
+		    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle) * timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle) * timeMult;
+
+		patchW[cellIndex] = StokesIFun :: W 
+		    (
+		        waveHeight_, 
+			RealwaterDepth_, 
+			waveKx, 
+			xGroup[cellGroup[cellIndex]-1], 
+			waveKy, 
+			yGroup[cellGroup[cellIndex]-1], 
+			waveOmega, 
+			currTime, 
+			wavePhase_, 
+			measuredLevels[cellGroup[cellIndex]-1]
+		);
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex]
+		    * timeMult;	
+	    }
+        }	
+    }

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesII.H

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = StokesIIFun :: U 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                waveKx, 
+                xGroup[cellGroup[cellIndex]-1], 
+                waveKy, 
+                yGroup[cellGroup[cellIndex]-1], 
+                waveOmega, 
+                currTime, 
+                wavePhase_, 
+                patchHeight[cellIndex]
+            );
+
+        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult;
+        patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult;
+
+        patchW[cellIndex] = StokesIIFun :: W 
+	    (
+                waveHeight_, 
+                RealwaterDepth_, 
+                waveKx, 
+                xGroup[cellGroup[cellIndex]-1], 
+                waveKy, 
+                yGroup[cellGroup[cellIndex]-1], 
+                waveOmega, 
+                currTime, 
+                wavePhase_, 
+                patchHeight[cellIndex]
+            );
+
+        patchW[cellIndex] = patchW[cellIndex]*timeMult;
+
+    }
+    else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = 0.0;
+        patchV[cellIndex] = 0.0;
+        patchW[cellIndex] = 0.0;
+    } 
+    else 
+    {	
+	if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    < measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		    - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = StokesIIFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle)*timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle)*timeMult;
+
+		patchW[cellIndex] = StokesIIFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	patchHeight[cellIndex]
+                    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+	    }
+	}
+	else if ( calculatedLevel[cellGroup[cellIndex]-1] 
+	    > measuredLevels[cellGroup[cellIndex]-1] )
+	{
+	    if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+	    {
+	        patchU[cellIndex] = StokesIIFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchV[cellIndex] = patchU[cellIndex] 
+		    * sin(waveAngle)*timeMult;
+		patchU[cellIndex] = patchU[cellIndex] 
+		    * cos(waveAngle)*timeMult;
+
+		patchW[cellIndex] = StokesIIFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                   );
+
+		patchW[cellIndex] = patchW[cellIndex]*timeMult;			
+	    }
+	    else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]))
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+		    - zInf[cellIndex];
+		auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+		auxiliarTotal = auxiliar/auxiliarTotal;
+
+		patchU[cellIndex] = StokesIIFun :: U 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * sin(waveAngle) * timeMult;
+		patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+		    * cos(waveAngle) * timeMult;
+
+		patchW[cellIndex] = StokesIIFun :: W 
+		    (
+                        waveHeight_, 
+                    	RealwaterDepth_, 
+                    	waveKx, 
+                    	xGroup[cellGroup[cellIndex]-1], 
+                    	waveKy, 
+                    	yGroup[cellGroup[cellIndex]-1], 
+                    	waveOmega, 
+                    	currTime, 
+                    	wavePhase_, 
+                    	measuredLevels[cellGroup[cellIndex]-1]
+                    );
+
+		patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+	    }
+        }		
+    }

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesV.H

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+    if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = stokesVFun :: U 
+	(
+            RealwaterDepth_, 
+            waveKx, 
+            waveKy, 
+            lambdaStokesV_, 
+            wavePeriod_, 
+            xGroup[cellGroup[cellIndex]-1], 
+            yGroup[cellGroup[cellIndex]-1], 
+            currTime, 
+            wavePhase_, 
+            patchHeight[cellIndex]
+        );
+
+        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult;
+        patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult;
+
+        patchW[cellIndex] = stokesVFun :: W 
+	(
+            RealwaterDepth_, 
+            waveKx, 
+            waveKy, 
+            lambdaStokesV_, 
+            wavePeriod_, 
+            xGroup[cellGroup[cellIndex]-1], 
+            yGroup[cellGroup[cellIndex]-1], 
+            currTime, 
+            wavePhase_, 
+            patchHeight[cellIndex]
+        );
+
+        patchW[cellIndex] = patchW[cellIndex]*timeMult;
+
+    }
+    else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1])
+    {
+        patchU[cellIndex] = 0.0;
+        patchV[cellIndex] = 0.0;
+        patchW[cellIndex] = 0.0;
+    } 
+    else 
+    {	
+        if ( calculatedLevel[cellGroup[cellIndex]-1] < 
+            measuredLevels[cellGroup[cellIndex]-1] )
+        {
+            if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+	            - zInf[cellIndex];
+	        auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+	        auxiliarTotal = auxiliar/auxiliarTotal;
+    
+	        patchU[cellIndex] = stokesVFun :: U 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+	            xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    patchHeight[cellIndex]
+                );
+
+	        patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+	            * sin(waveAngle)*timeMult;
+	        patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+	            * cos(waveAngle)*timeMult;
+
+	        patchW[cellIndex] = stokesVFun :: W 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+	    	    xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    patchHeight[cellIndex]
+                );
+
+		    patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		        * timeMult;	
+            }
+        }
+        else if ( calculatedLevel[cellGroup[cellIndex]-1] 
+            > measuredLevels[cellGroup[cellIndex]-1] )
+        {
+	    if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])
+            {
+	        patchU[cellIndex] = stokesVFun :: U 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_,    
+                    wavePeriod_, 
+	            xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    measuredLevels[cellGroup[cellIndex]-1]
+                );
+
+       	        patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) 
+		    * timeMult;
+                patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) 
+		    * timeMult;
+
+	        patchW[cellIndex] = stokesVFun :: W 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy,    
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+                    xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    measuredLevels[cellGroup[cellIndex]-1]
+                );
+
+	        patchW[cellIndex] = patchW[cellIndex]*timeMult;			
+	    }
+	    else if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) 
+	        & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) )
+	    {
+	        auxiliar = calculatedLevel[cellGroup[cellIndex]-1] 
+	            - zInf[cellIndex];
+	        auxiliarTotal = zSup[cellIndex]-zInf[cellIndex];
+	        auxiliarTotal = auxiliar/auxiliarTotal;
+
+	        patchU[cellIndex] = stokesVFun :: U 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+                    xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    measuredLevels[cellGroup[cellIndex]-1]
+                );
+    
+	        patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+	            * sin(waveAngle) * timeMult;
+	        patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] 
+	            * cos(waveAngle) * timeMult;
+
+	        patchW[cellIndex] = stokesVFun :: W 
+		(
+                    RealwaterDepth_, 
+                    waveKx, 
+                    waveKy, 
+                    lambdaStokesV_, 
+                    wavePeriod_, 
+                    xGroup[cellGroup[cellIndex]-1], 
+                    yGroup[cellGroup[cellIndex]-1], 
+                    currTime, 
+                    wavePhase_, 
+                    measuredLevels[cellGroup[cellIndex]-1]
+                );
+
+	        patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] 
+		    * timeMult;	
+	    }
+        }		
+    }

integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/velocityProfile.H

@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+   IH-Cantabria 2015 (http://www.ihcantabria.com/en/)
+   IHFOAM 2015 (http://ihfoam.ihcantabria.com/) 
+
+   Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
+               Gabriel Barajas   (barajasg@unican.es)
+\*---------------------------------------------------------------------------*/
+
+if ( waveType_ == "regular" )
+{
+    if ( waveTheory_ == "StokesI" )
+    {
+        #include "profileStokesI.H"
+        #include "profileGENAB.H"
+    }
+    else if ( waveTheory_ == "StokesII" )
+    {
+        #include "profileStokesII.H"
+        #include "profileGENAB.H"
+    }
+    else if ( waveTheory_ == "StokesV" )
+    {
+        #include "profileStokesV.H"
+        #include "profileGENAB.H"
+    }
+    else if ( waveTheory_ == "Cnoidal" )
+    {
+        #include "profileCnoidal.H"
+        #include "profileGENAB.H"
+    }
+}
+else if ( waveType_ == "solitary" )
+{
+    #include "profileBoussinesq.H"
+    #include "profileGENAB.H"
+}

integration/genAbs/waveGeneration/Make/files

@@ -0,0 +1,5 @@
+IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C
+IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C
+../common/waveFun.C
+
+LIB = $(FOAM_USER_LIBBIN)/libIHwaveGenerationV2esi

integration/genAbs/waveGeneration/Make/options

@@ -0,0 +1,16 @@
+sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
+sinclude $(RULES)/mplib$(WM_MPLIB)
+
+EXE_INC = \
+    -DOFVERSION=$(OF_VERSION) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I./IH_Waves_InletVelocity/velProfiles \
+    -I./IH_Waves_InletVelocityWind/velProfilesWind \
+    -I../common \
+    -I../common/checks \
+    -I../common/calculateWaterLevel \
+     $(PFLAGS) $(PINC)
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    $(PLIBS)

integration/genAbs/waveGeneration/localMake

@@ -0,0 +1,25 @@
+#!/bin/bash
+
+if [ $WM_PROJECT == "foam" ]; then
+    ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'`
+else
+    ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'`
+fi
+
+#IHC_dbg
+echo $ofversion
+#----
+
+export OF_VERSION=$ofversion
+
+wclean
+
+wmake libso
+
+if (( $? )) ; then 
+    echo "\n\nIH Wave generation boundary conditions (V2.0) compilation failed"
+    exit; else
+    echo -e "\n\nIH Wave generation boundary conditions (V2.0) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n";  
+fi
+
+wclean