/** \file PenaltyOperator.h */

#ifndef AMDIS_PENALTYOPERATOR_H
#define AMDIS_PENALTYOPERATOR_H

#include "ElementLevelSet.h"
#include "Flag.h"
#include "Operator.h"
#include "SurfaceOperator.h"
#include "ElementLevelSet.h"

namespace AMDiS {
    class ElInfo;
    class FiniteElemSpace;
}

using namespace AMDiS;
using namespace std;

// ===========================================================================
// ===== class PenaltyOperator ===============================================
// ===========================================================================
//
// Class Description:
//
// ===========================================================================
class PenaltyOperator : public Operator
{
public:
  /**
   * Constructor.
   */
  PenaltyOperator(Flag operatorType_,
		  ElementLevelSet *elLS_,
		  double factor_,
		  bool penaltyCoeffFlag_,
		  FiniteElemSpace *rowFESpace_,
		  FiniteElemSpace *colFESpace_ = NULL)
    : Operator(operatorType_, rowFESpace_, colFESpace_),
      elLS(elLS_),
      elStatus(ElementLevelSet::LEVEL_SET_UNDEFINED),
      factor(factor_),
      penaltyCoeffFlag(penaltyCoeffFlag_),
      surfaceOp(NULL),
      dim(getRowFESpace()->getMesh()->getDim()),
      degree(getRowFESpace()->getBasisFcts()->getDegree())
  {
    TEST_EXIT(elLS->getLevelSetFct() && elLS->getMesh())
      ("ElementLevelSet not initialized!\n");

    tempCoords = NEW VectorOfFixVecs<DimVec<double> >(dim, dim, NO_INIT);
  };

  /**
   * Destructor.
   */
  ~PenaltyOperator()
  {
    if (surfaceOp) {
      DELETE surfaceOp;
    }
/*     if (dim == 3) { */
      DELETE tempCoords;
/*     } */
  };

  /** 
   * Calculates the element matrix for this ElInfo and adds it multiplied by
   * factor to userMat.
   * In addition to \ref Operator::getElementMatrix(), it distinguishes
   * between elements divided by the boundary (level set zero intersects the 
   * element) and elements lying completely inside or outside the integration 
   * domain.
   */
  void getElementMatrix(const ElInfo *elInfo, 
			ElementMatrix& userMat, 
			double factor = 1.0);

  /** 
   * Calculates the element vector for this ElInfo and adds it multiplied by
   * factor to userVec.
   * Similar to \ref getElementMatrix it distinguishes between elements 
   * divided by the boundary and elements lying completely inside or outside 
   * the integration domain.
   */
  void getElementVector(const ElInfo *elInfo, 
			ElementVector& userVec, 
			double factor = 1.0);

protected:

  /**
   * Calculate coefficient of the penalty term.
   */
  double getPenaltyCoeff(const ElInfo *elInfo);

 protected:
  /**
   * Holds level set function and functionalities for intersection point
   * calculation.
   */
  ElementLevelSet *elLS;

  /**
   * Indicator for the position of an element. It shows whether an element 
   * lies inside the integration domain, outside of the integration domain 
   * or on the boundary.
   */
  int elStatus;

  /** 
   * Coefficient of the penalty term is 1/eps. eps depends on the size of
   * element and on factor.
   */
  double eps;

  /**
   * Factor needed for eps.
   */
  double factor;

  /**
   * Flag to indicate whether to use penalty coefficient.
   *   true - use penalty coefficient
   *   false - do not use
   */
  bool penaltyCoeffFlag;

  /**
   * Surface operator for surface integration
   */
  SurfaceOperator *surfaceOp;

  /**
   * Problem Dimension.
   */
  int dim;
 
  /**
   * Degree of basis functions.
   */
  int degree;
 
  /**
   * Variables used for calculation.
   */
  VectorOfFixVecs<DimVec<double> > *tempCoords;
};

#endif  // AMDIS_PENALTYOPERATOR_H