src/operator/OperatorOnUnknown.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program 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.
    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file OperatorOnUnknown.cpp
  \author E. Lunéville
  \since 4 mar 2012
  \date 9 may 2012

  \brief Implementation of xlifepp::OperatorOnUnknown class members and related functions
*/

//===============================================================================
// library dependencies
#include "OperatorOnUnknown.hpp"
//===============================================================================

namespace xlifepp
{

//------------------------------------------------------------------------------------
// OperatorOnUnknown member functions
//------------------------------------------------------------------------------------
// constructor of new OperatorOnUnknown
OperatorOnUnknown::OperatorOnUnknown(const Unknown* un, DiffOpType ty)
{
  difOp_p = findDifferentialOperator(ty);
  u_p = un;
  conjugateUnknown_ = false;
  if(un!=0 && un->conjugate())
    {
      conjugateUnknown_ = true;
      un->conjugate(false); // reset the conjugate state flag of unknown
    }
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = false;
  dimsRes_ = dimPair(0,0);
  setStructure();
}

OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, DiffOpType ty)
{
  difOp_p = findDifferentialOperator(ty);
  u_p = &un;
  conjugateUnknown_ = false;
  if(un.conjugate())
    {
      conjugateUnknown_ = true;
      un.conjugate(false); // reset the conjugate state flag of unknown
    }
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = false;
  dimsRes_ = dimPair(0,0);
  setStructure();
}

OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, DiffOpType ty, const std::vector<complex_t>& cs)
{
  difOp_p = findDifferentialOperator(ty);
  coefs_ = cs;
  u_p = &un;
  conjugateUnknown_ = false;
  if(un.conjugate())
    {
      conjugateUnknown_ = true;
      un.conjugate(false); // reset the conjugate state flag of unknown
    }
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = false;
  dimsRes_ = dimPair(0,0);
  setStructure();
}

// useful constructor : (fun op u) or (u op fun)
OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const Function& fun, AlgebraicOperator aop, bool lef)
{
  difOp_p = findDifferentialOperator(_id);
  u_p = &un;
  conjugateUnknown_ = false;
  if(un.conjugate())
    {
      conjugateUnknown_ = true;
      un.conjugate(false); // reset the conjugate state flag of unknown
    }
  dimsRes_ = dimPair(0,0);
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = lef;
  if(lef) { leftOperand_p = new Operand(fun, aop); }
  else { rightOperand_p = new Operand(fun, aop); }
  setStructure();
}

// useful constructor : (opfun op u) or (u op opfun)
OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const OperatorOnFunction& opfun, AlgebraicOperator aop, bool lef)
{
  difOp_p = findDifferentialOperator(_id);
  u_p = &un;
  conjugateUnknown_ = false;
  if(un.conjugate())
    {
      conjugateUnknown_ = true;
      un.conjugate(false); // reset the conjugate state flag of unknown
    }
  dimsRes_ = dimPair(0,0);
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = lef;
  if(lef) { leftOperand_p = new Operand(opfun, aop);}
  else { rightOperand_p = new Operand(opfun, aop);}
  setStructure();
}

// useful constructor : (val op u) or (u op val)
OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const Value& val, AlgebraicOperator aop, bool lef)
{
  difOp_p = findDifferentialOperator(_id);
  u_p = &un;
  conjugateUnknown_ = false;
  if(un.conjugate())
    {
      conjugateUnknown_ = true;
      un.conjugate(false); // reset the conjugate state flag of unknown
    }
  dimsRes_ = dimPair(0,0);
  leftOperand_p = 0;
  rightOperand_p = 0;
  leftPriority_ = lef;
  if(lef) { leftOperand_p = new Operand(val, aop); }
  else { rightOperand_p = new Operand(val, aop); }
  setStructure();
}

// copy constructor
OperatorOnUnknown::OperatorOnUnknown(const OperatorOnUnknown& opu)
{
  difOp_p = opu.difOp_p;
  conjugateUnknown_ = opu.conjugateUnknown_;
  u_p = opu.u_p;
  type_ = opu.type_;
  struct_ = opu.struct_;
  dimsRes_ = opu.dimsRes_;
  coefs_=opu.coefs_;
  leftOperand_p = 0;
  rightOperand_p = 0; // full copy
  if(opu.leftOperand_p != 0) { leftOperand_p = new Operand(*opu.leftOperand_p); }
  if(opu.rightOperand_p != 0) { rightOperand_p = new Operand(*opu.rightOperand_p); }
  leftPriority_ = opu.leftPriority_;
}

// assignment operator
OperatorOnUnknown& OperatorOnUnknown::operator=(const OperatorOnUnknown& opu)
{
  difOp_p = opu.difOp_p;
  conjugateUnknown_ = opu.conjugateUnknown_;
  u_p = opu.u_p;
  type_ = opu.type_;
  struct_ = opu.struct_;
  dimsRes_ = opu.dimsRes_;
  leftPriority_ = opu.leftPriority_;
  coefs_=opu.coefs_;
  if(leftOperand_p != 0) { delete leftOperand_p; }
  if(opu.leftOperand_p != 0) {leftOperand_p = new Operand(*opu.leftOperand_p);}
  else {leftOperand_p =0;}
  if(rightOperand_p != 0) { delete rightOperand_p; }
  if(opu.rightOperand_p != 0) {rightOperand_p = new Operand(*opu.rightOperand_p);}
  else {rightOperand_p =0;}
  return *this;
}

// destructor
OperatorOnUnknown::~OperatorOnUnknown()
{
  if(leftOperand_p != 0) { delete leftOperand_p; }
  if(rightOperand_p != 0) { delete rightOperand_p; }
}

// set returned type, structure, dims of returned value and check consistency
//   ### no longer used, replace by setStructure(), to be removed in next version ###
void OperatorOnUnknown::setReturnedType(const Unknown* un, DiffOpType ty)
{
  type_ = _real;
  struct_= _scalar;
  dimsRes_ = dimPair(1,1);
  if(un==0) return;

  type_=un->space()->valueType();   //may be complex in case of spectral space with complex functions

  dimen_t dimu = un->nbOfComponents();     // dimension of unknown
  dimen_t dimf = un->space()->dimFun();    // dimension of shape functions
  dimen_t dimv = std::max(dimu,dimf);      // dimension for vector case
  dimen_t dimd = un->space()->dimDomain(); // dimension of space domain
  dimen_t dims = un->space()->dimPoint();  // dimension of points, may be greater than dimd

  dimsRes_ = dimPair(dimu,1);

  switch(ty)
    {
      case _grad:
      case _gradS:
      case _gradG:
        {
          struct_ = _vector;
          dimsRes_=dimPair(dimd,1);
          if(dimv > 1) { struct_ = _matrix; dimsRes_ = dimPair(dimd, dimv); }
          break;
        }
      case _div:
      case _divS:
      case _divG:
        {
          struct_ = _scalar;
          dimsRes_ =dimPair(1,1);
          break;
        }
      case _curl:
      case _curlS:
      case _curlG:
        {
          struct_ = _vector;
          dimsRes_ = dimPair(dimv,1);
          break; // what about scalar curl in 2D
        }
      case _ntimes:
        {
          struct_ = _vector;
          dimsRes_ = dimPair(dimv,1);
          if(dimv > 1)
            {
              where("OperatorOnUnknown::setReturnedType(...)");
              error("operator_baddiffop", words("diffop", ty));
            }
          break;
        }
      case _ndot:
        {
          struct_ = _scalar;
          dimsRes_ = dimPair(1,1);
          if(dimv != dims)
            {
              where("OperatorOnUnknown::setReturnedType(...)");
              error("operator_baddiffop", words("diffop", ty));
            }
          break;
        }
      case _ncross:
      case _ncrossncross:
      case _ncrosscurl:
      case _ndiv:
        {
          struct_ = _vector;
          dimsRes_ = dimPair(dims,1);
          break;
        }
      case _ndotgrad:
        {
          struct_ = _scalar;
          dimsRes_ = dimPair(1,1);
          if(dimv > 1) { struct_ = _vector; dimsRes_ = dimPair(dimv,1);}
          break;
        }
      case _ncrossgrad:
        {
          struct_ = _scalar;
          dimsRes_ = dimPair(1,1);
          if(dims > 2) { struct_ = _vector; dimsRes_ = dimPair(dims,1);}
          break;
        }
      case _epsilon:
      case _epsilonG:
        {
          struct_ = _matrix;
          dimsRes_ = dimPair(dimv,dimv);
          break;
        }
      case _epsilonR:
        {
          struct_ = _vector;
          dimsRes_ = dimPair((dimv*(dimv+1))/2,1);
          break;
        }
      case _voigtToM:
        {
          struct_ = _matrix;
          dimsRes_ = dimPair(dimv,dimv);
          break;
        }
      default :
        {
          struct_ = _scalar;
          dimsRes_ = dimPair(1,1);
          if(dimv > 1) { struct_ = _vector; dimsRes_ = dimPair(dimv,1);}
        }
    }
}

// set returned type, structure, dims of returned value
void OperatorOnUnknown::setStructure()
{
  type_=_real;
  struct_=_scalar;
  dimsRes_=dimPair(0,0);
  if(u_p==0) return;   //no operator
  // set value type
  type_=u_p->space()->valueType();     //may be complex in case of spectral space with complex functions
  if(type_==_real && leftOperand_p!=0)  type_=leftOperand_p->valueType();
  if(type_==_real && rightOperand_p!=0) type_=rightOperand_p->valueType();
  //retry dimensions from unknown
  dimen_t dimu = u_p->nbOfComponents();     // dimension of unknown
  dimen_t dimf = u_p->space()->dimFun();    // dimension of shape functions
  dimen_t dimv = std::max(dimu,dimf);       // dimension for vector case
  dimen_t dims = u_p->space()->dimPoint();  // dimension of points, may be greater than dimd
  if(dimv > 1 && dimv < dims) dimv = dims;  // shapevalues will be mapped from dimv to dims
  Vector<real_t> np(1,1.);
  if(dims==2) {np=Vector<real_t>(2,1.);}
  else if(dims==3) {np=Vector<real_t>(3,1.);}
  Point P=np, Q=P+1.;
  dimen_t d=0, m=0;
  //evaluate operator on a fake shape function
  if(type_==_real)
    {
      std::vector<real_t> sv(dimv,1.);
      Vector<real_t> val;
      std::vector< std::vector<real_t> > dsv(std::max(dims,dimv),sv);
      if(!hasFunction())   eval(sv, dsv, dimv, val, d, m,&np);
      else if(hasFun())
        {
          if(normalRequired()) setNx(&np);
          eval(P, sv, dsv, dimv, val, d, m, &np);
        }
      else if(hasKernel())
        {
          if(xnormalRequiredByOpKernel()) setNx(&np);
          if(ynormalRequiredByOpKernel()) setNy(&np);
          eval(P, Q, sv, dsv, dimv, val, d, m, &np, &np);
        }
    }
  else // complex case
    {
      std::vector<real_t> sv(dimv,1.);
      std::vector< std::vector<real_t> > dsv(std::max(dims,dimv),sv);
      Vector<complex_t> val;
      if(!hasFunction())   eval(sv, dsv, dimv, val, d, m,&np);
      else if(hasFun())
        {
          if(normalRequired()) setNx(&np);
          eval(P, sv, dsv, dimv, val, d, m, &np);
        }
      else if(hasKernel())
        {
          if(xnormalRequiredByOpKernel()) setNx(&np);
          if(ynormalRequiredByOpKernel()) setNy(&np);
          eval(P, Q, sv, dsv, dimv, val, d, m, &np, &np);
        }
    }
  dimsRes_= dimPair(d/m,m);
  if(dimsRes_.first==1) struct_=_scalar;
  else if(dimsRes_.second==1) struct_=_vector;
  else struct_=_matrix;
  return;
}

//! order of polynom used by unknown (degree of basis)
number_t OperatorOnUnknown::unknownDegree() const
{
  if(unknown()->space()->typeOfSpace()==_feSpace)
    return unknown()->space()->feSpace()->maxDegree();  // new way, shape type no longer used
  return 5; //not a FE space (no polynoms) (default order is chosen as 5)
}

/*! order of polynom used by unknown (degree od basis), decreased by order of derivation
    as we use numtype (order of interpolation) we distinguish Pk or Qk using a ShapeType argument*/
number_t OperatorOnUnknown::degree() const
{
  if(unknown()->space()->typeOfSpace()==_feSpace)
      {
          number_t d = unknown()->space()->feSpace()->maxDegree(); // degree of shape functions
          if(hasFunction()) d*=2;  //degree is double to take into account the function or kernel involved
          d-=diffOrder();          //take into derivative
          if(d==0) d=1;            //min degree is 1
          return d;
      }
  return 5;   //not a FE space (no polynoms) (default order is chosen as 5)
}

//! order of differential operator involved
dimen_t OperatorOnUnknown::diffOrder() const
{
  if(difOp_p->type()==_epsilonG) return dimen_t(coefs_[0].real());   //special case
  return difOp_p->order();
}

//! true if operator involves some functions or kernels
bool OperatorOnUnknown::hasFunction() const
{
  if(leftOperand_p != 0 && (leftOperand_p->isFunction()  || leftOperand_p->isKernel())) { return true; }
  if(rightOperand_p != 0 && (rightOperand_p->isFunction() || rightOperand_p->isKernel())) { return true; }
  return false;
}

//! true if operator involves some functions or kernels
bool OperatorOnUnknown::hasKernel() const
{
  if(leftOperand_p != 0  &&  leftOperand_p->isKernel()) { return true; }
  if(rightOperand_p != 0 &&  rightOperand_p->isKernel()) { return true; }
  return false;
}

//! true if operator involves some functions
bool OperatorOnUnknown::hasFun() const
{
  if(leftOperand_p != 0  &&  leftOperand_p->isFunction()) { return true; }
  if(rightOperand_p != 0 &&  rightOperand_p->isFunction()) { return true; }
  return false;
}

//! direct acces to Kernel (unique Function, left priority)
const Function* OperatorOnUnknown::functionp() const
{
  if(leftOperand_p != 0  &&  leftOperand_p->isFunction()) { return leftOperand_p->functionp(); }
  if(rightOperand_p != 0 &&  rightOperand_p->isFunction()) { return rightOperand_p->functionp(); }
  return 0;
}

//! direct acces to Kernel (unique Kernel, left priority)
const Kernel* OperatorOnUnknown::kernelp() const
{
  if(leftOperand_p != 0  &&  leftOperand_p->isKernel()) { return leftOperand_p->kernelp(); }
  if(rightOperand_p != 0 &&  rightOperand_p->isKernel()) { return rightOperand_p->kernelp(); }
  return 0;
}

//! direct acces to Kernel (unique Kernel, left priority)
const OperatorOnKernel* OperatorOnUnknown::opkernelp() const
{
  if(leftOperand_p != 0  &&  leftOperand_p->isKernel()) { return leftOperand_p->opkernelp(); }
  if(rightOperand_p != 0 &&  rightOperand_p->isKernel()) { return rightOperand_p->opkernelp(); }
  return 0;
}


//! change Kernel in operator (unique Kernel, left priority)
void OperatorOnUnknown::changeKernel(Kernel* ker)
{
  if(leftOperand_p != 0  &&  leftOperand_p->isKernel()) {leftOperand_p->changeKernel(ker); return;}
  if(rightOperand_p != 0 &&  rightOperand_p->isKernel()) {rightOperand_p->changeKernel(ker); return;}
  where("OperatorOnUnknown::changeKernel");
  error("null_pointer","kernel");
}

//! true if operator involves left or right operand
bool OperatorOnUnknown::hasOperand() const
{
  return (leftOperand_p != 0 || rightOperand_p !=0);
}

//! true if normal is required in computation
bool OperatorOnUnknown::normalRequired() const
{
  if(difOp_p->normalRequired()) return true;
  if(leftOperand_p != 0 && leftOperand_p->normalRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->normalRequired()) return true;
  return false;
}

//! true if normal is required in computation
bool OperatorOnUnknown::xnormalRequired() const
{
  if(difOp_p->normalRequired()) return true;
  if(leftOperand_p != 0 && leftOperand_p->xnormalRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->xnormalRequired()) return true;
  return false;
}

//! true if normal is required in computation
bool OperatorOnUnknown::ynormalRequired() const
{
  if(difOp_p->normalRequired()) return true;
  if(leftOperand_p != 0 && leftOperand_p->ynormalRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->ynormalRequired()) return true;
  return false;
}

//! true if normal is required in computation by operator on Kernel
bool OperatorOnUnknown::xnormalRequiredByOpKernel() const
{
  if(leftOperand_p != 0 && leftOperand_p->xnormalRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->xnormalRequired()) return true;
  return false;
}

//! true if normal is required in computation by operator on Kernel
bool OperatorOnUnknown::ynormalRequiredByOpKernel() const
{
  if(leftOperand_p != 0 && leftOperand_p->ynormalRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->ynormalRequired()) return true;
  return false;
}

//!< true if GeomElemnt is required by operator
bool OperatorOnUnknown::elementRequired() const
{
  if(leftOperand_p != 0 && leftOperand_p->elementRequired()) return true;
  if(rightOperand_p != 0 && rightOperand_p->elementRequired()) return true;
  return false;
}

//! update value, structure and dimensions of result of (operator)(op)(valueType,strucType)
// dims : dimensions of value
void OperatorOnUnknown::updateReturnedType(AlgebraicOperator op, ValueType vt, StrucType st, dimPair dims, bool left)
{
  if(vt == _complex || type_ == _complex) { type_ = _complex; } //update value type

  switch(op)  // update structure type
    {
      case _product :
        {
          switch(struct_) // current opu struct
            {
              case _scalar :
                {
                  struct_ = st;  // scalar*any
                  return;
                }
              case _vector :
                {
                  struct_ = _vector;
                  if(st == _scalar)  // vector*scalar or scalar*vector
                    {
                      dimsRes_= dims;
                      return;
                    }
                  if(st == _matrix) // matrix*vector or vector*matrix
                    {
                      if(left) dimsRes_= dimPair(dims.first,1);
                      else dimsRes_= dimPair(dims.second,1);
                      return;
                    }
                  break;
                }
              case _matrix :
                {
                  struct_ = _matrix;
                  if(st == _scalar)  return; // scalar*matrix
                  if(st == _vector) // matrix*vector
                    {
                      if(left) dimsRes_= dimPair(dims.first,dimsRes_.second);
                      else dimsRes_= dimPair(dimsRes_.second, dims.second);
                      return;
                    }
                  if(st == _matrix) // matrix*matrix
                    {
                      if(left) dimsRes_= dimPair(dimsRes_.second,1);
                      else dimsRes_= dimPair(dimsRes_.first, 1);
                      return;
                    }
                  break;
                }
              default:
                break;
            }
          break;
        }
      case _innerProduct :
        {
          switch(struct_) // current opu struct
            {
              case _scalar :
                {
                  if(st == _scalar) // scalar.scalar
                    {
                      struct_ = _scalar;
                      return;
                    }
                  break;
                }
              case _vector :
                {
                  if(st == _vector) // vector.vector
                    {
                      struct_ = _scalar;
                      dimsRes_=dimPair(1,1);
                      return;
                    }
                  break;
                }
              case _matrix :
                break;
              default:
                break;
            }
          break;
        }
      case _crossProduct :
        {
          switch(struct_) // current opu struct
            {
              case _scalar :
              case _matrix :
                break;
              case _vector :
                {
                  if(st == _vector) // vector^vector
                    {
                      if(dims.first<=2) //return a scalar
                        {
                          struct_ = _scalar;
                          dimsRes_=dimPair(1,1);
                          return;
                        }
                      struct_ = _vector;
                      dimsRes_=dimPair(dims.first,1);
                      return;
                    }
                }
              default:
                break;
            }
        }
        break;
      case _contractedProduct :
        {
          switch(struct_) // current opu struct
            {
              case _scalar :
              case _vector :
                break;
              case _matrix :
                {
                  if(st == _matrix) // matrix%matrix
                    //if matrix have the same size return a scalar (standard contracted product)
                    //else generalized contracted product returning a matrix of same size of original (see computation)
                    {
                      if(dimsRes_==dims)
                        {
                          struct_ = _scalar;
                          dimsRes_=dimPair(1,1);
                          return;
                        }
                      struct_ = _matrix;
                      return;
                    }
                }
              default:
                break;
            }
        }
    }
  if(left)
    {
      where("OperatorOnUnknown::updateReturnedType(...)");
      error("operator_badop", words("structure", struct_), words("algop", op), words("structure", st));
    }
  else
    {
      where("OperatorOnUnknown::updateReturnedType(...)");
      error("operator_badop", words("structure", st), words("algop", op), words("structure", struct_));
    }
}

//------------------------------------------------------------------------------------
// external functions relates to OperatorOnUnknown
//------------------------------------------------------------------------------------

// user interface functions
// be cautious with memory, OperatorOnUnknown are created on dynamic memory
// at the end of the chaining construction, be sure to free it safely

//==========================================================================================
//  operators involving DifferentialOperator and Unknown (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown& id(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _id));
}
OperatorOnUnknown& d0(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d0));
}
OperatorOnUnknown& dt(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d0));
}
OperatorOnUnknown& d1(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d1));
}
OperatorOnUnknown& dx(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d1));
}
OperatorOnUnknown& d2(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d2));
}
OperatorOnUnknown& dy(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d2));
}
OperatorOnUnknown& d3(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d3));
}
OperatorOnUnknown& dz(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _d3));
}
OperatorOnUnknown& grad(const Unknown& un)
{
  return *(new OperatorOnUnknown(&un, _grad));
}
OperatorOnUnknown& nabla(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _grad));
}
OperatorOnUnknown& div(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _div));
}
OperatorOnUnknown& curl(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _curl));
}
OperatorOnUnknown& rot(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _curl));
}
OperatorOnUnknown& gradS(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _gradS));
}
OperatorOnUnknown& nablaS(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _gradS));
}
OperatorOnUnknown& divS(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _divS));
}
OperatorOnUnknown& curlS(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _curlS));
}
OperatorOnUnknown& rotS(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _curlS));
}
OperatorOnUnknown& nx(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ntimes));
}
OperatorOnUnknown& ndot(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ndot));
}
OperatorOnUnknown& ndotgrad(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ndotgrad));
}
OperatorOnUnknown& ncrossgrad(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ncrossgrad));
}
OperatorOnUnknown& ncross(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ncross));
}
OperatorOnUnknown& ncrossncross(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ncrossncross));
}
OperatorOnUnknown& ncrossntimes(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ncrossntimes));
}
//OperatorOnUnknown& ncrossndot(const Unknown& un)
//{
//  return*(new OperatorOnUnknown(&un, _ncrossndot));
//}
OperatorOnUnknown& ntimesndot(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ntimesndot));
}
OperatorOnUnknown& ndiv(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ndiv));
}
OperatorOnUnknown& ncrosscurl(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _ncrosscurl));
}
OperatorOnUnknown& epsilon(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _epsilon));
}

OperatorOnUnknown& epsilonR(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _epsilonR));
}

OperatorOnUnknown& voigtToM(const Unknown& un)   //special operator : convert Voigt vector to tensor matrix
{
  return*(new OperatorOnUnknown(&un, _voigtToM));
}

// generalized differential operator
OperatorOnUnknown& gradG(const Unknown& un, const complex_t& ax, const complex_t& ay,
                         const complex_t& az, const complex_t& at)
{
  std::vector<complex_t> cs(4);
  cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
  OperatorOnUnknown* opu = new OperatorOnUnknown(un, _gradG, cs);
  return *opu;
}

OperatorOnUnknown& nablaG(const Unknown& un, const complex_t& ax, const complex_t& ay,
                          const complex_t& az, const complex_t& at)
{
  return gradG(un, ax, ay, az, at);
}

OperatorOnUnknown& divG(const Unknown& un, const complex_t& ax, const complex_t& ay,
                        const complex_t& az, const complex_t& at)
{
  std::vector<complex_t> cs(4);
  cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
  OperatorOnUnknown* opu = new OperatorOnUnknown(un, _divG, cs);
  return *opu;
}

OperatorOnUnknown& curlG(const Unknown& un, const complex_t& ax, const complex_t& ay,
                         const complex_t& az, const complex_t& at)
{
  std::vector<complex_t> cs(4);
  cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
  OperatorOnUnknown* opu = new OperatorOnUnknown(un, _curlG, cs);
  return *opu;
}

OperatorOnUnknown& rotG(const Unknown& un, const complex_t& ax, const complex_t& ay,
                        const complex_t& az, const complex_t& at)
{
  return curlG(un, ax, ay, az, at);
}

OperatorOnUnknown& epsilonG(const Unknown& un, const complex_t& id, const complex_t& ax,
                            const complex_t& ay, const complex_t& az)
{
  std::vector<complex_t> cs(4);
  cs[0]=id; cs[1]=ax; cs[2]=ay; cs[3]=az;
  OperatorOnUnknown* opu = new OperatorOnUnknown(un, _epsilonG, cs);
  return *opu;
}

OperatorOnUnknown& mean(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _mean));
}

OperatorOnUnknown& jump(const Unknown& un)
{
  return*(new OperatorOnUnknown(&un, _jump));
}

//==========================================================================================
// operators involving UnitaryVector and Unknown (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown& operator*(UnitaryVector v, const Unknown& un)       // n*u same as nx(u)
{
  if(v == _n) return *(new OperatorOnUnknown(&un, _ntimes));
  if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossntimes));
  where("UnitaryVector * Unknown");
  error("operator_unexpected", words("diffop",v)+ "? * u");
  return *(new OperatorOnUnknown()); // dummy return
}

OperatorOnUnknown& operator|(UnitaryVector v, const Unknown& un)       // n|u same as ndot(u)
{
  if(v == _n) return *(new OperatorOnUnknown(&un, _ndot));
//  if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossndot));
  where("UnitaryVector | Unknown");
  error("operator_unexpected", words("diffop",v)+ "? | u");
  return *(new OperatorOnUnknown()); // dummy return
}

OperatorOnUnknown& operator^(UnitaryVector v, const Unknown& un)       // n^u same as ncross(u)
{
  if(v == _n) return *(new OperatorOnUnknown(&un, _ncross));
  if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossncross));
  where("UnitaryVector ^ Unknown");
  error("operator_unexpected", words("diffop",v)+ "? ^ u");
  return *(new OperatorOnUnknown()); // dummy return
}

OperatorOnUnknown& operator^(const Unknown& un, UnitaryVector v)       // u^n same as -ncross(u)
{
  OperatorOnUnknown* opu=0;
  if(v == _n)       opu = new OperatorOnUnknown(&un, _ncross);
  if(v == _ncrossn) opu = new OperatorOnUnknown(&un, _ncrossncross);
  where("Unknown ^ UnitaryVector");
  if(opu==0) error("operator_unexpected", words("diffop",v)+ "u ^ ?");
  opu->leftOperand() = new Operand(-1.,_product);
  return *opu;
}

OperatorOnUnknown& operator*(const Unknown& un, UnitaryVector v)       // u*n same as nx(u)
{return v * un;}

OperatorOnUnknown& operator|(const Unknown& un, UnitaryVector v)       // u|n same as ndot(u)
{return v | un;}

OperatorOnUnknown& operator|(UnitaryVector v, OperatorOnUnknown& opu) // n|grad(u) same as ndotgrad(u)
{
  if(v == _n && opu.difOpType() == _grad) {opu.difOp_() = findDifferentialOperator(_ndotgrad); return opu;}
  if(v == _n && opu.difOpType() == _id)   {opu.difOp_() = findDifferentialOperator(_ndot); return opu;}
//  if(v == _ncrossn && opu.difOpType() == _id)   {opu.difOp_() = findDifferentialOperator(_ncrossndot); return opu;}
  where("UnitaryVector | OperatorOnUnknown");
  error("operator_unexpected", words("diffop",v)+ " | "+ words("diffop",opu.difOpType()));
  return opu;
}

OperatorOnUnknown& operator*(UnitaryVector v, OperatorOnUnknown& opu) // n*, n*div, n^n*, n*n|
{
  if(v == _n && opu.difOpType() == _id)  {opu.difOp_() = findDifferentialOperator(_ntimes); return opu;}
  if(v == _n && opu.difOpType() == _div) {opu.difOp_() = findDifferentialOperator(_ndiv); return opu;}
  if(v == _ncrossn && opu.difOpType() == _id)   {opu.difOp_() = findDifferentialOperator(_ncrossntimes); return opu;}
  if(v == _n && opu.difOpType()==_ndot) {opu.difOp_() = findDifferentialOperator(_ntimesndot); return opu;}
  where("UnitaryVector * OperatorOnUnknown");
  error("operator_unexpected", words("diffop",v)+ " * "+ words("diffop",opu.difOpType()));
  return opu;
}

OperatorOnUnknown& operator^(UnitaryVector v, OperatorOnUnknown& opu) // n^(n^u) or n^curl(u) or n^grad(u) same as ncrossncross(u) or ncrosscurl(u) or ncrossgrad(u)
{
  DiffOpType ty = opu.difOpType();
  if(v != _n)
    {
      where("UnitaryVector ^ OperatorOnUnknown");
      error("operator_unexpected", words("normal",v)+ "? ^ "+ words("diffop",ty));
    }
  switch(ty)
    {
      case _ncross: opu.difOp_() = findDifferentialOperator(_ncrossncross); break;
      case _grad: opu.difOp_() = findDifferentialOperator(_ncrossgrad); break;
      case _curl: opu.difOp_() = findDifferentialOperator(_ncrosscurl); break;
      default:
        where("UnitaryVector ^ OperatorOnUnknown");
        error("operator_unexpected", words("normal",v)+ " ^ "+ words("diffop",ty)+"?");
    }
  return opu;
}

OperatorOnUnknown& operator^(OperatorOnUnknown& opu, UnitaryVector v) // (n^u)^n or curl(u)^n or grad(u)^n same as -ncrossncross(u) or -ncrosscurl(u) or -ncrossgrad(u)
{
  DiffOpType ty = opu.difOpType();
  if(v != _n)
    {
      where("OperatorOnUnknown ^ UnitaryVector");
      error("operator_unexpected", words("diffop",ty)+ " ^ " + words("normal",v) + "?");
    }
  switch(ty)
    {
      case _ncross: opu.difOp_() = findDifferentialOperator(_ncrossncross); break;
      case _grad: opu.difOp_() = findDifferentialOperator(_ncrossgrad); break;
      case _curl: opu.difOp_() = findDifferentialOperator(_ncrosscurl); break;
      default:
        where("OperatorOnUnknown ^ UnitaryVector");
        error("operator_unexpected", words("diffop",ty)+ " ^ "+words("normal",v)+ "?");
    }
  if(opu.leftOperand()==0) opu.leftOperand()= new Operand(-1.,_product);
  else
    {
      if(opu.leftOperand()->isValue()) opu.leftOperand()->value().opposite();
      else
        {
          where("OperatorOnUnknown ^ UnitaryVector");
          error("opposite_operand_impossible");
        }
    }
  return opu;
}

OperatorOnUnknown& operator|(OperatorOnUnknown& opu, UnitaryVector v) // only grad(u)|n same as ndotgrad(u)
{return v | opu;}

OperatorOnUnknown& operator*(OperatorOnUnknown& opu, UnitaryVector v) // only div*n same as ndiv(u)
{return v * opu;}

//==========================================================================================
// operators involving Unknown and Function (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown&  operator*(const Unknown& un, const Function& fun)      // u*F
{return *(new OperatorOnUnknown(un, fun, _product, false));}

OperatorOnUnknown&  operator|(const Unknown& un, const Function& fun)      // u|F
{return *(new OperatorOnUnknown(un, fun, _innerProduct, false));}

OperatorOnUnknown&  operator^(const Unknown& un, const Function& fun)      // u^F
{return *(new OperatorOnUnknown(un, fun, _crossProduct, false));}

OperatorOnUnknown&  operator%(const Unknown& un, const Function& fun)       // u%F
{return *(new OperatorOnUnknown(un, fun, _contractedProduct, false));}

OperatorOnUnknown&  operator*(const Function& fun, const Unknown& un)       // F*u
{return *(new OperatorOnUnknown(un, fun, _product, true));}

OperatorOnUnknown&  operator|(const Function& fun, const Unknown& un)       // F|u
{return *(new OperatorOnUnknown(un, fun, _innerProduct, true));}

OperatorOnUnknown&  operator^(const Function& fun, const Unknown& un)       // F^u
{return *(new OperatorOnUnknown(un, fun, _crossProduct, true));}

OperatorOnUnknown& operator%(const Function& fun, const Unknown& un)        // F%u
{return *(new OperatorOnUnknown(un, fun, _contractedProduct, true));}

//==========================================================================================
// operators involving Unknown and OperatorOnFunction (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown&  operator*(const Unknown& un, const OperatorOnFunction& opfun)      // u*op(F)
{return *(new OperatorOnUnknown(un, opfun, _product, false));}

OperatorOnUnknown&  operator|(const Unknown& un, const OperatorOnFunction& opfun)      // u|op(F)
{return *(new OperatorOnUnknown(un, opfun, _innerProduct, false));}

OperatorOnUnknown&  operator^(const Unknown& un, const OperatorOnFunction& opfun)      // u^op(F)
{return *(new OperatorOnUnknown(un, opfun, _crossProduct, false));}

OperatorOnUnknown&  operator%(const Unknown& un, const OperatorOnFunction& opfun)       // u%op(F)
{return *(new OperatorOnUnknown(un, opfun, _contractedProduct, false));}

OperatorOnUnknown&  operator*(const OperatorOnFunction& opfun, const Unknown& un)       // op(F)*u
{return *(new OperatorOnUnknown(un, opfun, _product, true));}

OperatorOnUnknown&  operator|(const OperatorOnFunction& opfun, const Unknown& un)       // op(F)|u
{return *(new OperatorOnUnknown(un, opfun, _innerProduct, true));}

OperatorOnUnknown&  operator^(const OperatorOnFunction& opfun, const Unknown& un)       // op(F)^u
{return *(new OperatorOnUnknown(un, opfun, _crossProduct, true));}

OperatorOnUnknown& operator%(const OperatorOnFunction& opfun, const Unknown& un)        // op(F)%u
{return *(new OperatorOnUnknown(un, opfun, _contractedProduct, true));}

//==========================================================================================
//  operators involving Unknown and Value (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown& operator*(const Unknown& un, const Value& val)          // u*val
{return *(new OperatorOnUnknown(un, val, _product, false));}

OperatorOnUnknown& operator|(const Unknown& un, const Value& val)          // u|val
{return *(new OperatorOnUnknown(un, val, _innerProduct, false));}

OperatorOnUnknown& operator^(const Unknown& un, const Value& val)          // u^val
{return *(new OperatorOnUnknown(un, val, _crossProduct, false));}

OperatorOnUnknown& operator%(const Unknown& un, const Value& val)          // u%val
{return *(new OperatorOnUnknown(un, val, _contractedProduct, false));}

OperatorOnUnknown& operator*(const Value& val, const Unknown& un)          // val*u
{return *(new OperatorOnUnknown(un, val, _product, true));}

OperatorOnUnknown& operator|(const Value& val, const Unknown& un)          // val|u
{return *(new OperatorOnUnknown(un, val, _innerProduct, true));}

OperatorOnUnknown& operator^(const Value& val, const Unknown& un)          // val^u
{return *(new OperatorOnUnknown(un, val, _crossProduct, true));}

OperatorOnUnknown& operator%(const Value& val, const Unknown& un)          // val%u
{return *(new OperatorOnUnknown(un, val, _contractedProduct, true));}

//==========================================================================================
//  operators involving Unknown and explicit value (create new OperatorOnUnknown)
//==========================================================================================
OperatorOnUnknown& operator*(const Unknown& un, const real_t& val)
{return fromUnknownVal(un,val,_product);}

OperatorOnUnknown& operator|(const Unknown& un, const real_t& val)
{return fromUnknownVal(un,val,_innerProduct);}

OperatorOnUnknown& operator^(const Unknown& un, const real_t& val)
{return fromUnknownVal(un,val,_crossProduct);}

OperatorOnUnknown& operator%(const Unknown& un, const real_t& val)
{return fromUnknownVal(un,val,_contractedProduct);}

OperatorOnUnknown& operator*(const real_t& val, const Unknown& un)
{return fromValUnknown(un,val,_product);}

OperatorOnUnknown& operator|(const real_t& val, const Unknown& un)
{return fromValUnknown(un,val,_innerProduct);}

OperatorOnUnknown& operator^(const real_t& val, const Unknown& un)
{return fromValUnknown(un,val,_crossProduct);}

OperatorOnUnknown& operator%(const real_t& val, const Unknown& un)
{return fromValUnknown(un,val,_contractedProduct);}

OperatorOnUnknown& operator*(const Unknown& un, const complex_t& val)
{return fromUnknownVal(un,val,_product);}

OperatorOnUnknown& operator|(const Unknown& un, const complex_t& val)
{return fromUnknownVal(un,val,_innerProduct);}

OperatorOnUnknown& operator^(const Unknown& un, const complex_t& val)
{return fromUnknownVal(un,val,_crossProduct);}

OperatorOnUnknown& operator%(const Unknown& un, const complex_t& val)
{return fromUnknownVal(un,val,_contractedProduct);}

OperatorOnUnknown& operator*(const complex_t& val, const Unknown& un)
{return fromValUnknown(un,val,_product);}

OperatorOnUnknown& operator|(const complex_t& val, const Unknown& un)
{return fromValUnknown(un,val,_innerProduct);}

OperatorOnUnknown& operator^(const complex_t& val, const Unknown& un)
{return fromValUnknown(un,val,_crossProduct);}

OperatorOnUnknown& operator%(const complex_t& val, const Unknown& un)
{return fromValUnknown(un,val,_contractedProduct);}


//==========================================================================================
// operators involving OperatorOnUnknown and Function (update current OperatorOnUnknown)
//==========================================================================================
// update left and right operation with Function :  op(u) aop Function or Function aop op(u)
OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const Function& fun, AlgebraicOperator aop)
{
  if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
  //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), false);
  opu.rightOperand() = new Operand(fun, aop);
  if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
  opu.setStructure();
  return opu;
}
OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const Function& fun, AlgebraicOperator aop)
{
  if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
  //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), true);
  opu.leftOperand() = new Operand(fun, aop);
  if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
  opu.setStructure();
  return opu;
}

// left and right operation with Function :  op(u) aop F or F aop op(u)
OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Function& fun) // d(u) * F
{
  return updateRight(opu, fun, _product);
}
OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Function& fun) // d(u) | F
{
  return updateRight(opu, fun, _innerProduct);
}
OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Function& fun) // d(u) ^ F
{
  return updateRight(opu, fun, _crossProduct);
}
OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Function& fun) // d(u) % F
{
  return updateRight(opu, fun, _contractedProduct);
}
OperatorOnUnknown& operator*(const Function& fun, OperatorOnUnknown& opu) // F * d(u)
{
  return updateLeft(opu, fun, _product);
}
OperatorOnUnknown& operator|(const Function& fun, OperatorOnUnknown& opu) // F | d(u)
{
  return updateLeft(opu, fun, _innerProduct);
}
OperatorOnUnknown& operator^(const Function& fun, OperatorOnUnknown& opu) // F ^ d(u)
{
  return updateLeft(opu, fun, _crossProduct);
}
OperatorOnUnknown& operator%(const Function& fun, OperatorOnUnknown& opu) // F % d(u)
{
  return updateLeft(opu, fun, _contractedProduct);
}

//==========================================================================================
// operators involving OperatorOnUnknown and OperatorOnFunction (update current OperatorOnUnknown)
//==========================================================================================
// update left and right operation with Function :  op(u) aop op(F) or op(F) aop op(u)
OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const OperatorOnFunction& fun, AlgebraicOperator aop)
{
  if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
  //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), false);
  opu.rightOperand() = new Operand(fun, aop);
  if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
  opu.setStructure();
  return opu;
}
OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const OperatorOnFunction& fun, AlgebraicOperator aop)
{
  if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
  //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), true);
  opu.leftOperand() = new Operand(fun, aop);
  if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
  opu.setStructure();
  return opu;
}

// left and right operation with OperatorOnFunction :  op(u) aop F or F aop op(u)
OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) * F
{
  return updateRight(opu, fun, _product);
}
OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) | F
{
  return updateRight(opu, fun, _innerProduct);
}
OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) ^ F
{
  return updateRight(opu, fun, _crossProduct);
}
OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) % F
{
  return updateRight(opu, fun, _contractedProduct);
}
OperatorOnUnknown& operator*(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F * d(u)
{
  return updateLeft(opu, fun, _product);
}
OperatorOnUnknown& operator|(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F | d(u)
{
  return updateLeft(opu, fun, _innerProduct);
}
OperatorOnUnknown& operator^(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F ^ d(u)
{
  return updateLeft(opu, fun, _crossProduct);
}
OperatorOnUnknown& operator%(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F % d(u)
{
  return updateLeft(opu, fun, _contractedProduct);
}

//==========================================================================================
// operators involving OperatorOnUnknown and scalar (update current OperatorOnUnknown)
//==========================================================================================
// left and right operation with real scalar :  op(u) aop r or r aop op(u)
OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const real_t& val)
{return updateRight(opu, val, _product);}

OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const real_t& val)
{return updateRight(opu, val, _innerProduct);}

OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const real_t& val)
{return updateRight(opu, val, _crossProduct);}

OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const real_t& val)
{return updateRight(opu, val, _contractedProduct);}

OperatorOnUnknown& operator*(const real_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _product);}

OperatorOnUnknown& operator|(const real_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _innerProduct);}

OperatorOnUnknown& operator^(const real_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _crossProduct);}

OperatorOnUnknown& operator%(const real_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _contractedProduct);}

// left and right operation with complex scalar :  op(u) aop c or c aop op(u)
OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const complex_t& val)
{return updateRight(opu, val, _product);}

OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const complex_t& val)
{return updateRight(opu, val, _innerProduct);}

OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const complex_t& val)
{return updateRight(opu, val, _crossProduct);}

OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const complex_t& val)
{return updateRight(opu, val, _contractedProduct);}

OperatorOnUnknown& operator*(const complex_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _product);}

OperatorOnUnknown& operator|(const complex_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _innerProduct);}

OperatorOnUnknown& operator^(const complex_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _crossProduct);}

OperatorOnUnknown& operator%(const complex_t& val, OperatorOnUnknown& opu)
{return updateLeft(opu, val, _contractedProduct);}

//==========================================================================================
// operators involving OperatorOnUnknown and Value (update current OperatorOnUnknown)
//==========================================================================================
// update left and right operation with Value :  op(u) aop Value or Value aop op(u)
OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const Value& val, AlgebraicOperator o)
{
  if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
  //opu.updateReturnedType(o, val.valueType(), val.strucType(), val.dims(), false);
  opu.rightOperand() = new Operand(&val, o);
  if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
  opu.setStructure();
  return opu;
}
OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const Value& val, AlgebraicOperator o)
{
  if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
  //opu.updateReturnedType(o, val.valueType(), val.strucType(), val.dims(), true);
  opu.leftOperand() = new Operand(&val, o);
  if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
  opu.setStructure();
  return opu;
}
OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Value& val) // d(u) * val
{
  return updateRight(opu, val, _product);
}
OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Value& val) // d(u) | val
{
  return updateRight(opu, val, _innerProduct);
}
OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Value& val) // d(u) ^ val
{
  return updateRight(opu, val, _crossProduct);
}
OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Value& val) // d(u) % val
{
  return updateRight(opu, val, _contractedProduct);
}
OperatorOnUnknown& operator*(const Value& val, OperatorOnUnknown& opu) // val * d(u)
{
  return updateLeft(opu, val, _product);
}
OperatorOnUnknown& operator|(const Value& val, OperatorOnUnknown& opu) // val | d(u)
{
  return updateLeft(opu, val, _innerProduct);
}
OperatorOnUnknown& operator^(const Value& val, OperatorOnUnknown& opu) // val ^ d(u)
{
  return updateLeft(opu, val, _crossProduct);
}
OperatorOnUnknown& operator%(const Value& val, OperatorOnUnknown& opu) // val % d(u)
{
  return updateLeft(opu, val, _contractedProduct);
}

//==========================================================================================
// various facilities
//==========================================================================================

// check is operation between two operatoronunknown is consistant in terms of structure
// the result have to be a scalar
bool checkConsistancy(const OperatorOnUnknown& opu, AlgebraicOperator aop, const OperatorOnUnknown& opv)
{
  if(opu.strucType() == _scalar && opv.strucType() == _scalar) { return true; }
  switch(aop)
    {
      case _product : // only scalar * scalar
        if(opu.strucType() == _scalar && opv.strucType() == _scalar) { return true; }
        break;
      case _innerProduct : // only vector | vector
        if(opu.strucType() == _vector && opv.strucType() == _vector) { return true; }
        break;
      case _contractedProduct : // only matrix % matrix
        if(opu.strucType() == _matrix && opv.strucType() == _matrix) { return true; }
        break;
      case _crossProduct : // never
        break;
    }
  return false;
}


// compare OperatorOnUnknown (same unknowns, same diff operators, same functions ...)
bool operator==(const OperatorOnUnknown& opu, const OperatorOnUnknown& opv)
{
  if(opu.unknown()!=opv.unknown()) return false;
  if(opu.conjugateUnknown()!=opv.conjugateUnknown()) return false;
  if(opu.difOp()!=opv.difOp()) return false;
  if(opu.leftOperand()!=0 && opv.leftOperand()!=0)
    if(*opu.leftOperand()!=*opv.leftOperand()) return false;
  if(opu.rightOperand()!=0 && opv.rightOperand()!=0)
    if(*opu.rightOperand()!=*opv.rightOperand()) return false;
  if(opu.coefs()!=opu.coefs()) return false;
  return true;
}

bool operator!=(const OperatorOnUnknown& opu, const OperatorOnUnknown& opv)
{
  return !(opu==opv);
}

//==========================================================================================
// printing facilities
//==========================================================================================
void OperatorOnUnknown::print(std::ostream& os) const
{
  if(theVerboseLevel > 0)
    {
      os << words("operator") << " " << difOp_p->name() << " "<<u_p->name() ;
      if(conjugateUnknown_) { os << "(" << words("conjugate") << ")"; }
      os << " " << words("returns a") << " " << words("value", type_) << " " << words("structure", struct_);
      if(struct_!=_scalar)   os << " ("<<dimsRes_.first<<" x "<<dimsRes_.second<<")";
      if(coefs_.size()>0) os<<"  coefficients : "<<coefs_;
      if(leftOperand_p != 0)  { os << "\n   " << words("left operand") << " : " << *leftOperand_p; }
      if(rightOperand_p != 0) { os << "\n   " << words("right operand") << " : " << *rightOperand_p; }
      if(leftOperand_p != 0 && rightOperand_p != 0)
        {
          if(leftPriority_) { os << " (" << words("left priority") << ")"; }
          else              { os << " (" << words("right priority") << ")"; }
        }
      os << "\n";
    }
}

string_t OperatorOnUnknown::asString() const
{
  string_t s="";
  if(leftOperand_p != 0)     s += leftOperand_p->asString();
  if(difOp_p->type()!=_id)   s += difOp_p->name()+"("+u_p->name()+")";
  else                       s += u_p->name();
  if(rightOperand_p != 0)    s += rightOperand_p->asString();
  return s;
}

void OperatorOnUnknown::printsymbolic(std::ostream& os) const
{
  if(leftOperand_p != 0)  { leftOperand_p->printsymbolic(os);}
  if(difOp_p->type()!=_id)
    os<<difOp_p->name()<<"("<<u_p->name()<<")";
  else
    os<<u_p->name();
  if(rightOperand_p != 0)  { rightOperand_p->printsymbolic(os);}
}

std::ostream& operator<<(std::ostream& os, const OperatorOnUnknown& opu)
{
  opu.print(os);
  return os;
}

} // end of namespace xlifepp