/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

Description
    A pair of surfaces in contact.

\*---------------------------------------------------------------------------*/

#include "contactPatchPair.H"
#include "contactProblem.H"


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

// Construct from dictionary
Foam::contactPatchPair::contactPatchPair
(
    const word& name,
    contactProblem& cp,
    const dictionary& dict
)
:
    contactActive_(dict.lookup("contactActive")),
    name_(name),
    cp_(cp),
    alg_(intersection::algorithmNames_.read(dict.lookup("projectionAlgo"))),
    dir_(intersection::directionNames_.read(dict.lookup("projectionDir"))),
    masterPatch_(dict.lookup("masterPatch"), cp.mesh().boundaryMesh()),
    slavePatch_(dict.lookup("slavePatch"), cp.mesh().boundaryMesh()),
    contactGapTol_(readScalar(dict.lookup("contactGapTol"))),
    penaltyScaleFactor_(readScalar(dict.lookup("penaltyScaleFactor"))),
    returnScaleFactor_(readScalar(dict.lookup("returnScaleFactor"))),
    contactIterNum_(0),
    contactFilePtr_(NULL),
    slaveFaceZoneName_(slavePatch_.name() + "FaceZone"),
    masterFaceZoneName_(masterPatch_.name() + "FaceZone"),
    slavePointZoneName_(slavePatch_.name() + "PointZone"),
    totalSlavePointForce_(
  			  cp_.mesh().boundaryMesh()[slavePatch_.index()].pointNormals().size(), vector::zero),
			  // IOobject
			  // (
			  //  "totalSlavePointForce",
			  //  cp_.mesh().time().timeName(),
			  //  cp_.mesh(),
			  //  IOobject::READ_IF_PRESENT,
			  //  IOobject::AUTO_WRITE
			  //  ),
    			  // vectorField(cp_.mesh().boundaryMesh()[slavePatch_.index()].pointNormals().size(), vector::zero)
			  // ),
    slavePointPenetration_(
			   cp_.mesh().boundaryMesh()[slavePatch_.index()].pointNormals().size(),
			   0.0
			   ),
    masterPointPenetration_(
			    cp_.mesh().boundaryMesh()[masterPatch_.index()].pointNormals().size(),
			    0.0
			    ),
    globalSlavePointPenetration_(
				 cp_.mesh().pointZones()[cp_.mesh().faceZones().findZoneID(slaveFaceZoneName_)].size(),
				 0.0
				 ),
    globalMasterPointPenetration_(
				  cp_.mesh().pointZones()[cp_.mesh().faceZones().findZoneID(masterFaceZoneName_)].size(),
				  0.0
				  ),
    oldTotalSlavePointForce_(
    			     cp_.mesh().boundaryMesh()[slavePatch_.index()].pointNormals().size(),
    			     vector::zero
    			     ),
    oldMinSlavePointPenetration_(0.0),
    penetrationLimit_(readScalar(dict.lookup("penetrationLimit"))),
    rigidMaster_(dict.lookup("rigidMaster")),
    interpolationMethod_(dict.lookup("interpolationMethod")),
    faceZoneMasterToSlaveInterpolator_(
				       cp_.mesh().faceZones()[cp_.mesh().faceZones().findZoneID(masterFaceZoneName_)](), // from
				       cp_.mesh().faceZones()[cp_.mesh().faceZones().findZoneID(slaveFaceZoneName_)](), // to zone
				       alg_,
				       dir_
				       ),
    faceZoneSlaveToMasterInterpolatorPtr_(NULL),
    faceZoneGgiInterpolatorPtr_(NULL),
    slaveInterpolator_(cp_.mesh().boundaryMesh()[slavePatch_.index()])
{
  Info << "\tConstructing contact patch pair"
       << "\n\t\tmaster patch:\t" << masterPatch_.name()
       << "\n\t\tslave patch:\t" << slavePatch_.name()
       << endl;
  
  if(rigidMaster_)
    Info << "\t\tThe master surface is considered rigid and is set as traction-free" << endl;

  //- Check interpolation scheme for passing tractions from slave to master
  if(interpolationMethod_ == "patchToPatch")
    {
      Info << "\t\tMethod for interpolation of traction from slave to master: patchToPatch" << endl;
      label masterFaceZoneID = cp_.mesh().faceZones().findZoneID(masterFaceZoneName_);
      label slaveFaceZoneID = cp_.mesh().faceZones().findZoneID(slaveFaceZoneName_);
      faceZoneSlaveToMasterInterpolatorPtr_ =
	new zoneToZoneInterpolation
	(
	 cp_.mesh().faceZones()[slaveFaceZoneID](), // from zone
	 cp_.mesh().faceZones()[masterFaceZoneID](), // to zone
	 alg_,
	 dir_
	 );
    }
  else if(interpolationMethod_ == "ggi")
    {
      Info << "\t\tMethod for interpolation of traction from slave to master: ggi" << endl;
      label masterFaceZoneID = cp_.mesh().faceZones().findZoneID(masterFaceZoneName_);
      label slaveFaceZoneID = cp_.mesh().faceZones().findZoneID(slaveFaceZoneName_);
      faceZoneGgiInterpolatorPtr_ =
	new ggiZoneInterpolation
	(
	 cp_.mesh().faceZones()[masterFaceZoneID](), // master zone
	 cp_.mesh().faceZones()[slaveFaceZoneID](), // slave zone
	 tensorField(0),
	 tensorField(0),
	 vectorField(0),
	 0.0,
	 0.0,
	 true,
	 ggiInterpolation::AABB
	 );
    }
  else
    {
      SeriousError << "\n\nTraction interpolation method '"
		   << interpolationMethod_ << "' not found!\n"
		   << "\nValid interpolation methods are:\n"
		   << "ggi\n"
		   << "patchToPatch"
		   << endl
		   << exit(FatalError);
    }

  //- only the master should create a contactInfo file
  if(Pstream::master())
    {
      contactFilePtr_ = new OFstream(name_);
      OFstream& contactFile = *contactFilePtr_;
      int width = 20;
      contactFile << "Time";
      contactFile.width(width);
      contactFile << "ContactIter";
      contactFile.width(width);
      contactFile << "PenaltyScale";
      contactFile.width(width);
      contactFile << "slaveContactVer";
      contactFile.width(width);
      contactFile << "penetration";
      contactFile.width(width);
      contactFile << "slaveVerReduce";
      contactFile.width(width);
      contactFile << "maxSlaveTrac" << endl;
    }
  Info << "\tContact patch pair constructed" << endl;
}
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //



// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //


Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::masterTouchFraction() const
{
  //--DOES NOT WORK IN PARALLEL YET BUT IS ONLY FOR VISUALISATION--/

  const fvMesh& mesh = cp_.mesh();

  scalarField pDistToInter(cp_.mesh().boundaryMesh()[masterPatch_.index()].size(), 0.0);
  //    = slaveToMasterInterpolateDeformed.pointDistanceToIntersection();
  
  scalarField vertexMasterGap =  pDistToInter;
 
// Calculate area in contact
    const faceList& masterPatchLocalFaces =
        mesh.boundaryMesh()[masterPatch_.index()].localFaces();

    const pointField& masterPatchLocalPoints =
        mesh.boundaryMesh()[masterPatch_.index()].localPoints();

    tmp<scalarField> ttouchFrac
    (
        new scalarField(masterPatchLocalFaces.size(), 0)
    );
    scalarField& touchFrac = ttouchFrac();

    forAll (masterPatchLocalFaces, faceI)
    {
        touchFrac[faceI] =
            masterPatchLocalFaces[faceI].areaInContact
            (
                masterPatchLocalPoints,
                vertexMasterGap
            );
    }

    return ttouchFrac;
}




Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::slaveTouchFraction() const
{
  //--DOES NOT WORK IN PARALLEL YET BUT IS ONLY FOR VISUALISATION--/

  const fvMesh& mesh = cp_.mesh();

  scalarField vertexSlaveGap = slavePointPenetration_;
  //    (cp_.mesh().boundaryMesh()[slavePatch_.index()].size(), 0.0);
    
  // Calculate area in contact
    const faceList& slavePatchLocalFaces =
      mesh.boundaryMesh()[slavePatch_.index()].localFaces();

    const pointField& slavePatchLocalPoints =
        mesh.boundaryMesh()[slavePatch_.index()].localPoints();

    tmp<scalarField> ttouchFrac
    (
        new scalarField(slavePatchLocalFaces.size(), 0)
    );
    scalarField& touchFrac = ttouchFrac();

    forAll (slavePatchLocalFaces, faceI)
    {
        touchFrac[faceI] =
            slavePatchLocalFaces[faceI].areaInContact
            (
                slavePatchLocalPoints,
                vertexSlaveGap
            );
    }

    return ttouchFrac;
}




Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::masterGapPoints() const
{
  scalarField globalMasterPointPenetration
    (
     cp_.mesh().boundaryMesh()[masterPatch_.index()].meshPoints().size(),
     0.0
     );

  scalarField vertexMasterGap = masterPointPenetration_;
  //globalMasterPointPenetration;
   
    tmp<scalarField> tcontactGapPoints
    (
	  new scalarField(vertexMasterGap.size(), 0)
    );
    scalarField& contactGapPoints = tcontactGapPoints();

    forAll (vertexMasterGap, pointI)
    {
      contactGapPoints[pointI] = vertexMasterGap[pointI];
    }

    return tcontactGapPoints;
}




Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::slaveGapPoints() const
{
  const scalarField& slavePointPenetration = slavePointPenetration_;

  tmp<scalarField> tcontactGapPoints
    (
     new scalarField(slavePointPenetration.size(), 0)
    );
    scalarField& contactGapPoints = tcontactGapPoints();

    forAll (slavePointPenetration, pointI)
    {
      contactGapPoints[pointI] = slavePointPenetration[pointI];
    }

    return tcontactGapPoints;
}




Foam::tmp<Foam::vectorField>
Foam::contactPatchPair::masterPointForce() const
{
  //- returns zero at the moment
  vectorField masterPointForce
    (
     cp_.mesh().boundaryMesh()[masterPatch_.index()].meshPoints().size(),
     vector::zero
     );
  
  tmp<vectorField> tcontactPointForce
    (
     new vectorField(masterPointForce.size(), vector::zero)
     );
  vectorField& contactPointForce = tcontactPointForce();
  
  forAll (contactPointForce, pointI)
    {
      contactPointForce[pointI] = masterPointForce[pointI];
    }
  
  return tcontactPointForce;
}




Foam::tmp<Foam::vectorField>
Foam::contactPatchPair::slavePointForce() const
{ 
  vectorField slavePointForce = totalSlavePointForce_;
  
  tmp<vectorField> tcontactPointForce
    (
     new vectorField(slavePointForce.size(), vector::zero)
     );
  
  vectorField& contactPointForce = tcontactPointForce();

  forAll (contactPointForce, pointI)
    {
      contactPointForce[pointI] = slavePointForce[pointI];
    }
  
  return tcontactPointForce;
}




Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::masterContactPressure() const
{
    const fvMesh& mesh = cp_.mesh();

    vectorField n = mesh.boundary()[masterPatch_.index()].nf();

    const volSymmTensorField& sigma =
        mesh.objectRegistry::lookupObject<volSymmTensorField>("sigma");

    scalarField masterPressure
      = n & (n & sigma.boundaryField()[masterPatch_.index()]);

    tmp<scalarField> tcontactPressure
    (
        new scalarField(masterPressure.size(), 0)
    );
    scalarField& contactPressure = tcontactPressure();

    forAll (masterPressure, faceI)
    {
      contactPressure[faceI] = max(-1*masterPressure[faceI],0);
    }

    return tcontactPressure;
}




Foam::tmp<Foam::scalarField>
Foam::contactPatchPair::slaveContactPressure() const
{
    const fvMesh& mesh = cp_.mesh();

    vectorField n = mesh.boundary()[slavePatch_.index()].nf();

    const volSymmTensorField& sigma =
        mesh.objectRegistry::lookupObject<volSymmTensorField>("sigma");

    scalarField slavePressure
      = n & (n & sigma.boundaryField()[slavePatch_.index()]);

    tmp<scalarField> tcontactPressure
    (
        new scalarField(slavePressure.size(), 0)
    );
    scalarField& contactPressure = tcontactPressure();

    forAll (slavePressure, faceI)
    {
      contactPressure[faceI] = max(-1*slavePressure[faceI],0);
    }

    return tcontactPressure;
}

// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //


// ************************************************************************* //