xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/term/SuTermVector.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 SuTermVector.cpp
  \authors D. Martin, E. Lunéville
  \since 03 apr 2012
  \date  13 jan 2014

  \brief Implementation of xlifepp::SuTermVector class members and related utilities
 */

#include "SuTermVector.hpp"
#include "BilinearFormAsLinearForm.hpp"
#include "TermMatrix.hpp"
#include "geometry.h"
#include "utils.h"
#include <numeric>


namespace xlifepp
{
//===============================================================================
//member functions of SuTermVector class
//===============================================================================
//constructors, destructor, assignment
SuTermVector::SuTermVector(SuLinearForm* sulf, const string_t& na, bool noass)
{
  termType_ = _sutermVector;
  sulf_p = sulf;
  name_ = na;
  if(sulf_p != 0)  u_p = sulf_p->unknown();
  else u_p = 0;
  space_p = 0;
  computingInfo_.noAssembly = noass;
  entries_p = 0;
  scalar_entries_p = 0;
  if(sulf_p != 0) buildSubspaces();   //construct the required subspaces
}

SuTermVector::SuTermVector(const string_t& na, const Unknown* u, Space* sp,
                           ValueType vt, number_t n, dimen_t nv, bool noass)
{
  termType_ = _sutermVector;
  name_ = na;
  sulf_p = 0;
  u_p = u;
  space_p = sp;
  computingInfo_.noAssembly = noass;
  entries_p = 0;
  scalar_entries_p = 0;
  if(n > 0) entries_p = new VectorEntry(vt, nv, n);
}

//copy constructor (full copy of entries)
SuTermVector::SuTermVector(const SuTermVector& suV)
{
  copy(suV);
}

//full copy of entries assuming SutermVector is clean
void SuTermVector::copy(const SuTermVector& suV)
{
  termType_ = _sutermVector;
  sulf_p = suV.sulf_p;
  name_ = suV.name_;
  u_p = suV.u_p;
  space_p = suV.space_p;
  subspaces = suV.subspaces;
  computingInfo_ = suV.computingInfo_;
  entries_p = 0;
  scalar_entries_p = 0;
  if(suV.entries_p != 0) entries_p = new VectorEntry(*suV.entries_p);
  if(suV.scalar_entries_p != suV.entries_p && suV.scalar_entries_p != 0)
    {
      scalar_entries_p = new VectorEntry(*suV.scalar_entries_p);
    }
  cdofs_ = suV.cdofs_;
}

//SuTermVector constructor from Unknown and function (no linear form). Specialization of Function or OperatorOnFunction
//compute f(Mi) for all Mi,Lagrange Dof of space
template <>
SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, const OperatorOnFunction& opf, const string_t& na, bool noass)
{
  initFromFunction(u, dom, opf, na, noass);
}

template <>
SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, const Function& f, const string_t& na, bool noass)
{
  initFromFunction(u, dom, OperatorOnFunction(f), na, noass);
}

SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, funSR_t& f, const string_t& na, bool noass)
{
  Function F(f);
  OperatorOnFunction opf(F);
  initFromFunction(u, dom, opf, na, noass);
}

SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, funSC_t& f, const string_t& na, bool noass)
{
  Function F(f);
  OperatorOnFunction opf(F);
  initFromFunction(u, dom, opf, na, noass);
}

SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, funVR_t& f, const string_t& na, bool noass)
{
  Function F(f);
  OperatorOnFunction opf(F);
  initFromFunction(u, dom, opf, na, noass);
}

SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, funVC_t& f, const string_t& na, bool noass)
{
  Function F(f);
  OperatorOnFunction opf(F);
  initFromFunction(u, dom, opf, na, noass);
}

//specialization of SuTermVector constructor from Unknown and SymbolicFunction
template <>
SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, const SymbolicFunction& fs, const string_t& na, bool noass)
{
  if(u.nbOfComponents() > 1)
    {
      where("SuTermVector::SuTermVector(Unknown, DOmain, SymbolicFunction, String, bool)");
      error("unknown_not_vector");
    }
  termType_ = _sutermVector;
  computingInfo_.noAssembly = noass;
  name_ = na;
  sulf_p = 0;
  u_p = &u;
  //find subspace or create it
  Space* spu = u.space();
  space_p = spu->findSubSpace(&dom, spu);
  if(space_p == 0) space_p = new Space(dom, *spu, spu->name() + "_" + dom.name());
  //allocate entries
  number_t n = space_p->nbDofs();
  entries_p = 0;
  scalar_entries_p = 0;
  if(n == 0) return;   // no entries
  //fill entries with fs values
  ValueType vt = fs.valueType();
  entries_p = new VectorEntry(vt, _scalar, n, 1);
  if(vt == _real)
    {
      Vector<real_t>::iterator ite = entries_p->rEntries_p->begin();
      for(number_t i = 1; i <= n; i++, ite++) *ite = fs(space_p->dof(i).coords());
    }
  else
    {
      Vector<complex_t>::iterator ite = entries_p->cEntries_p->begin();
      for(number_t i = 1; i <= n; i++, ite++) *ite = fs(space_p->dof(i).coords());
    }
  computingInfo_.isComputed = true;
}

//specialization of SuTermVector constructor from Unknown and VariableName
template <>
SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, const VariableName& vn, const string_t& na, bool noass)
{
  if(u.nbOfComponents() > 1)
    {
      where("SuTermVector::SuTermVector(Unknown, Domain, VariableName, String, bool)");
      error("unknown_not_vector");
    }
  termType_ = _sutermVector;
  computingInfo_.noAssembly = noass;
  name_ = na;
  sulf_p = 0;
  u_p = &u;
  //find subspace or create it
  Space* spu = u.space();
  space_p = spu->findSubSpace(&dom, spu);
  if(space_p == 0) space_p = new Space(dom, *spu, spu->name() + "_" + dom.name());
  //allocate entries
  number_t n = space_p->nbDofs();
  entries_p = 0;
  scalar_entries_p = 0;
  if(n == 0) return;   // no entries
  //create symbolic function
  SymbolicFunction fs(vn);
  //fill entries with fs values
  ValueType vt = fs.valueType();
  entries_p = new VectorEntry(vt, _scalar, n, 1);
  if(vt == _real)
    {
      Vector<real_t>::iterator ite = entries_p->rEntries_p->begin();
      for(number_t i = 1; i <= n; i++, ite++) *ite = fs(space_p->dof(i).coords());
    }
  else
    {
      Vector<complex_t>::iterator ite = entries_p->cEntries_p->begin();
      for(number_t i = 1; i <= n; i++, ite++) *ite = fs(space_p->dof(i).coords());
    }
  computingInfo_.isComputed = true;
}

/*! apply function to SutermVector
    build new SutermVector Y from SuTermVector X with Y(i)=f(X(i))
    function f has to be consistent with SuTermVector coefficients
        use a scalar function for scalar SuTermVector and vector function for vector SuTermVector
*/
SuTermVector::SuTermVector(const SuTermVector& X, funSR1_t& f, const string_t& na)
{
  if(X.strucType() != _scalar || X.valueType() != _real)
    {
      where("SuTermVector::SuTermVector(SuTermVector, funSR1_t, String)");
      error("value_wrongtype", words("structure", X.strucType()), words("value", X.valueType()),
            words("structure", _scalar), words("value", _real));
    }
  copy(X);
  name_ = na;
  if(entries_p == 0 || entries_p->rEntries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, funSR1_t, String)");
      error("term_no_entries");
    }
  Vector<real_t>::iterator itv = entries_p->rEntries_p->begin();
  for(; itv != entries_p->rEntries_p->end(); ++itv) *itv = f(*itv);
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

SuTermVector::SuTermVector(const SuTermVector& X, funSC1_t& f, const string_t& na)
{
  if(X.strucType() != _scalar || X.valueType() != _complex)
    {
      where("SuTermVector::SuTermVector(SuTermVector, funSC1_t, String)");
      error("value_wrongtype", words("structure", X.strucType()), words("value", X.valueType()),
            words("structure", _scalar), words("value", _complex));
    }
  copy(X);
  name_ = na;
  if(entries_p == 0 || entries_p->cEntries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, funSC1_t, String)");
      error("term_no_Entries");
    }
  Vector<complex_t>::iterator itv = entries_p->cEntries_p->begin();
  for(; itv != entries_p->cEntries_p->end(); ++itv) *itv = f(*itv);
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

/*! apply function to two SutermVector's
    build new SutermVector Z from SuTermVectors X and Y with Z(i)=f(X(i),Y(i))
*/
SuTermVector::SuTermVector(const SuTermVector& X, const SuTermVector& Y, funSR2_t& f, const string_t& na)
{
  if(X.space_p != Y.space_p)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSR2_t, String);");
      error("term_mismatch_spaces", X.space_p->name(), Y.space_p->name());
    }
  if(X.strucType() != _scalar || X.valueType() != _real)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSR2_t, String);");
      error("value_wrongtype", words("structure", X.strucType()), words("value", X.valueType()),
            words("structure", _scalar), words("value", _real));
    }
  if(X.entries_p == 0 || Y.entries_p == 0 || X.entries_p->rEntries_p == 0 || Y.entries_p->rEntries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSR2_t, String);");
      error("term_no_entries");
    }
  if(X.entries_p->rEntries_p->size() != Y.entries_p->rEntries_p->size())
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSR2_t, String);");
      error("entry_mismatch_dims", X.entries_p->rEntries_p->size(), Y.entries_p->rEntries_p->size());
    }
  copy(X);
  name_ = na;
  Vector<real_t>::iterator itv = entries_p->rEntries_p->begin(), itX = X.entries_p->rEntries_p->begin(), itY = Y.entries_p->rEntries_p->begin();
  for(; itv != entries_p->rEntries_p->end(); ++itv, ++itX, ++itY) *itv = f(*itX, *itY);
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

SuTermVector::SuTermVector(const SuTermVector& X, const SuTermVector& Y, funSC2_t& f, const string_t& na)
{
  if(X.space_p != Y.space_p)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSC2_t, String);");
      error("term_mismatch_spaces", X.space_p->name(), Y.space_p->name());
    }
  if(X.strucType() != _scalar || X.valueType() != _complex)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSC2_t, String)");
      error("value_wrongtype", words("structure", X.strucType()), words("value", X.valueType()),
            words("structure", _scalar), words("value", _complex));
    }
  if(X.entries_p == 0 || Y.entries_p == 0 || X.entries_p->cEntries_p == 0 || Y.entries_p->cEntries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSC2_t, String);");
      error("term_no_entries");
    }
  if(X.entries_p->cEntries_p->size() != Y.entries_p->cEntries_p->size())
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, funSC2_t, String);");
      error("entry_mismatch_dims", X.entries_p->cEntries_p->size(), Y.entries_p->cEntries_p->size());
    }
  copy(X);
  name_ = na;
  Vector<complex_t>::iterator itv = entries_p->cEntries_p->begin(), itX = X.entries_p->cEntries_p->begin(), itY = Y.entries_p->cEntries_p->begin();
  for(; itv != entries_p->cEntries_p->end(); ++itv, ++itX, ++itY) *itv = f(*itX, *itY);
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

//constructor from symbolic function of TermVector's
/*!
   build the SuTermVector R = fs(X) where fs is a symbolic function of x
   when X is a scalar SuTermVector, performs
         R(i)=fs(X(i))  i=1,n   result is a scalar SuTermVector with X unknown
   when X is a vector SuTermVector (dim=d), performs component by component
         R_k(i)=fs(X_k(i))  i=1,n, k=1,d, result is a vector SuTermVector with X unknown
*/
SuTermVector::SuTermVector(const SuTermVector& X, const SymbolicFunction& fs, const string_t& na)
{
  copy(X);
  name_ = na;
  bool isComplex = (fs.valueType()==_complex || X.valueType()==_complex);
  if(isComplex && X.valueType()==_real) toComplex();
  if(entries_p->rEntries_p != 0) //real scalar case
    {
      Vector<real_t>::iterator itx = X.entries_p->rEntries_p->begin();
      if(!isComplex)
        {
          Vector<real_t>::iterator itr = entries_p->rEntries_p->begin();
          for(; itr != entries_p->rEntries_p->end(); ++itx, ++itr) *itr = fs(*itx);
        }
      else
        {
          Vector<complex_t>::iterator itr = entries_p->cEntries_p->begin();
          for(; itr != entries_p->cEntries_p->end(); ++itx, ++itr) *itr = fs(complex_t(*itx));
        }
    }
  else if(entries_p->cEntries_p != 0) //complex scalar case
    {
      Vector<complex_t>::iterator itx = X.entries_p->cEntries_p->begin();
      Vector<complex_t>::iterator itr = entries_p->cEntries_p->begin();
      for(; itr != entries_p->cEntries_p->end(); ++itx,++itr) *itr = fs(*itx);
    }

  else  // real/complex vector case
    {
      Unknown* us = const_cast<Unknown*>(X.up());
      dimen_t nbc=us->nbOfComponents();
      for(dimen_t k=0; k<nbc; ++k)
        {
          SuTermVector sut = SuTermVector(X(*(*us)[k+1].asComponent()),fs); //compute kth component
          if(!isComplex)
            {
              Vector<Vector<real_t> >::iterator itr = entries_p->rvEntries_p->begin();
              Vector<real_t>::iterator itxk = sut.entries_p->rEntries_p->begin();
              for(; itr!=entries_p->rvEntries_p->end(); ++itr,++itxk)(*itr)[k]=*itxk;
            }
          else
            {
              Vector<Vector<complex_t> >::iterator itr = entries_p->cvEntries_p->begin();
              Vector<complex_t>::iterator itxk = sut.entries_p->cEntries_p->begin();
              for(; itr!=entries_p->cvEntries_p->end(); ++itr,++itxk)(*itr)[k]=*itxk;
            }
        }
    }
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

/*!
   build the SuTermVector R = fs(X1,X2) where fs is a symbolic function of x1_ x2_
   when X1 and X2 are scalar SuTermVector, performs
         R(i)=fs(X1(i),X2(i))  i=1,n   result is a scalar SuTermVector with X1 unknown
   when X1 is a vector SuTermVector (dim=d) and X2 a scalar SuTermvector, performs component by component
         R_k(i)=fs(X1_k(i),X2(i))  i=1,n, k=1,d, result is a vector SuTermVector with X1 unknown
   when X1 is a scalar SuTermVector and X2 a vector SuTermvector (dim=d), performs component by component
         R_k(i)=fs(X1(i),X2_k(i))  i=1,n, k=1,d, result is a vector SuTermVector with X2 unknown
   when X1 and X2 are both vector SuTermvector, computation is not available

   NB : unknown can be change later using changeUnknown function
*/
SuTermVector::SuTermVector(const SuTermVector& X1, const SuTermVector& X2, const SymbolicFunction& fs, const string_t& na)
{
  if(X1.space_p != X2.space_p)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SymbolicFunction, String);");
      error("term_mismatch_spaces", X1.space_p->name(), X2.space_p->name());
    }
  if(X1.entries_p == 0 || X2.entries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_no_entries");
    }

  bool fsIsComplex = (fs.valueType()==_complex);
  name_ = na;
  if(X1.strucType() == _scalar && X2.strucType() ==_scalar) //scalar case
    {

      if(X1.entries_p->rEntries_p != 0)  //X1 is real
        {
          Vector<real_t>::const_iterator itv1 = X1.entries_p->rEntries_p->begin();
          copy(X1);
          if(X2.entries_p->rEntries_p != 0)    //X2 is real
            {
              if(!fsIsComplex) //real result
                {
                  Vector<real_t>::iterator itv2 = X2.entries_p->rEntries_p->begin();
                  Vector<real_t>::iterator itvr = entries_p->rEntries_p->begin();
                  for(; itvr != entries_p->rEntries_p->end(); ++itv1, ++itv2, ++itvr) *itvr = fs(*itv1, *itv2);
                }
              else //complex result
                {
                  toComplex(); //move to complex
                  Vector<real_t>::iterator itv2 = entries_p->rEntries_p->begin();
                  Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
                  for(; itvr != entries_p->cEntries_p->end(); ++itv1, ++itv2, ++itvr) *itvr = fs(complex_t(*itv1), complex_t(*itv2));
                }
            }
          else if(X2.entries_p->cEntries_p != 0) //X2 is complex
            {
              toComplex(); //move to complex
              Vector<complex_t>::iterator itv2 = X2.entries_p->cEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itv1, ++itv2, ++itvr) *itvr = fs(complex_t(*itv1), *itv2);
            }
        }
      else // X1 is complex, result is complex
        {
          copy(X1);
          Vector<complex_t>::iterator itv1 = entries_p->cEntries_p->begin();
          if(X2.entries_p->rEntries_p != 0)  // X2 is real
            {
              Vector<real_t>::const_iterator itv2 = X2.entries_p->rEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itv1, ++itv2, ++itvr) *itv1 = fs(*itv1, complex_t(*itv2));
            }
          else
            {
              Vector<complex_t>::const_iterator itv2 = X2.entries_p->cEntries_p->begin();
              for(; itv1 != entries_p->cEntries_p->end(); ++itv1, ++itv2) *itv1 = fs(*itv1, *itv2);
            }
        }
      //remove scalar representation
      if(scalar_entries_p != 0)
        {
          delete scalar_entries_p;
          scalar_entries_p = 0;
        }
      return;  //end of scalar case
    }

  //vector case
  if(X1.strucType() != _scalar && X2.strucType() !=_scalar)  //two vector SuTermVector not available
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_unhandled_structures", words("structure", X1.strucType()), words("structure", X2.strucType()),fs.asString());
    }
  Unknown* us =0;
  dimen_t nbc=0;
  bool isComplex = false;
  bool X1isVector=true;
  if(X1.strucType() != _scalar) // X1 vector, X2 scalar
    {
      us = const_cast<Unknown*>(X1.up());
      nbc=us->nbOfComponents();
      copy(X1);
      isComplex = (fs.valueType()==_complex || X2.valueType()==_complex || X1.valueType()==_complex);
      if(isComplex && X1.valueType()==_real) toComplex();
    }
  else // X2 vector, X1 scalar
    {
      X1isVector=false;
      us = const_cast<Unknown*>(X2.up());
      nbc=us->nbOfComponents();
      copy(X2);
      isComplex = (fs.valueType()==_complex || X2.valueType()==_complex || X1.valueType()==_complex);
      if(isComplex && X2.valueType()==_real) toComplex();
    }
  for(dimen_t k=0; k<nbc; ++k)
    {
      SuTermVector sut;
      if(X1isVector) sut = SuTermVector(X1(*(*us)[k+1].asComponent()),X2,fs); //compute kth component
      else sut = SuTermVector(X1,X2(*(*us)[k+1].asComponent()),fs);           //compute kth component
      if(!isComplex)
        {
          Vector<Vector<real_t> >::iterator itvr = entries_p->rvEntries_p->begin();
          Vector<real_t>::iterator itvk = sut.entries_p->rEntries_p->begin();
          for(; itvr!=entries_p->rvEntries_p->end(); ++itvr,++itvk)(*itvr)[k]=*itvk;
        }
      else
        {
          Vector<Vector<complex_t> >::iterator itvr = entries_p->cvEntries_p->begin();
          Vector<complex_t>::iterator itvk = sut.entries_p->cEntries_p->begin();
          for(; itvr!=entries_p->cvEntries_p->end(); ++itvr,++itvk)(*itvr)[k]=*itvk;
        }
    }
  //remove scalar representation
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

/*!
   build the SuTermVector R = fs(X1,X2,X3) where fs is a symbolic function of x_1 x_2 x_3 :
     R(i)=fs(X1(i),X2(i),X3(i))  i=1,n
   X1, X2, X3 have to be three scalar SutermVectors defined on the same space
   the unknown of the SuTermVector result is the unknown of X1
   it can be change later using changeUnknown function
*/
SuTermVector::SuTermVector(const SuTermVector& X1, const SuTermVector& X2, const SuTermVector& X3,
                           const SymbolicFunction& fs, const string_t& na)
{
  if(X1.space_p != X2.space_p)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String);");
      error("term_mismatch_spaces", X1.space_p->name(), X2.space_p->name());
    }
  if(X1.space_p != X3.space_p)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String);");
      error("term_mismatch_spaces", X1.space_p->name(), X2.space_p->name());
    }
  if(X1.strucType() != _scalar)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_mismatch_structures", words("structure", X1.strucType()), words("structure", _scalar));
    }
  if(X2.strucType() != _scalar)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_mismatch_structures", words("structure", X2.strucType()), words("structure", _scalar));
    }
  if(X3.strucType() != _scalar)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_mismatch_structures", words("structure", X3.strucType()), words("structure", _scalar));
    }
  if(X1.entries_p == 0 || X2.entries_p == 0 || X3.entries_p == 0)
    {
      where("SuTermVector::SuTermVector(SuTermVector, SuTermVector, SuTermVector, SymbolicFunction, String)");
      error("term_no_entries");
    }

  name_ = na;
  bool fsIsComplex = (fs.valueType()==_complex);
  ValueType vt1 = X1.valueType(), vt2 = X2.valueType(), vt3 = X3.valueType();
  copy(X1);
  if((vt2==_complex || vt3==_complex || fsIsComplex) && vt1==_real) toComplex();

  if(vt1 == _real)
    {
      Vector<real_t>::const_iterator itv1 = X1.entries_p->rEntries_p->begin();
      if(vt2 == _real)
        {
          Vector<real_t>::const_iterator itv2 = X2.entries_p->rEntries_p->begin();
          if(vt3 == _real)
            {
              Vector<real_t>::const_iterator itv3 = X3.entries_p->rEntries_p->begin();
              if(!fsIsComplex)
                {
                  Vector<real_t>::iterator itvr = entries_p->rEntries_p->begin();
                  for(; itvr != entries_p->rEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(*itv1, *itv2, *itv3);
                }
              else
                {
                  Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
                  for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(complex_t(*itv1), complex_t(*itv2), complex_t(*itv3));
                }
            }
          else // X3 is complex
            {
              Vector<complex_t>::const_iterator itv3 = entries_p->cEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(complex_t(*itv1), complex_t(*itv2), *itv3);
            }
        }
      else // X2 is complex
        {
          Vector<complex_t>::const_iterator itv2 = entries_p->cEntries_p->begin();
          if(vt3 == _real)   //X3 is real
            {
              Vector<real_t>::const_iterator itv3 = X3.entries_p->rEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(complex_t(*itv1), *itv2, complex_t(*itv3));
            }
          else //X3 is complex
            {
              Vector<complex_t>::const_iterator itv3 = X3.entries_p->cEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(complex_t(*itv1), *itv2, *itv3);
            }
        }
    }
  else //X1 is complex
    {
      Vector<complex_t>::iterator itv1 = entries_p->cEntries_p->begin();
      if(vt2 == _real)  //X2 is real
        {
          Vector<real_t>::const_iterator itv2 = X2.entries_p->rEntries_p->begin();
          if(vt3 == _real)  //X3 is real
            {
              Vector<real_t>::const_iterator itv3 = X3.entries_p->rEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(*itv1, complex_t(*itv2), complex_t(*itv3));
            }
          else //X3 is complex
            {
              Vector<complex_t>::const_iterator itv3 = X3.entries_p->cEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(*itv1, complex_t(*itv2), *itv3);
            }
        }
      else //X2 is complex
        {
          Vector<complex_t>::const_iterator itv2 = X2.entries_p->cEntries_p->begin();
          if(vt3 == _real)  //X3 is real
            {
              Vector<real_t>::const_iterator itv3 = X3.entries_p->rEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(*itv1, *itv2, complex_t(*itv3));
            }
          else //X3 is complex
            {
              Vector<complex_t>::const_iterator itv3 = X3.entries_p->cEntries_p->begin();
              Vector<complex_t>::iterator itvr = entries_p->cEntries_p->begin();
              for(; itvr != entries_p->cEntries_p->end(); ++itvr, ++itv1, ++itv2, ++itv3) *itvr = fs(*itv1, *itv2, *itv3);
            }
        }
    }
  if(scalar_entries_p != 0)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
}

/*! construct vector SutermVector from 2 scalar SuTermVector's
    scalar SuTermVector's have to be defined on the same space
    vu has to be a vector unknown at least d=2 */
SuTermVector::SuTermVector(const Unknown& vu, const SuTermVector& X1, const SuTermVector& X2, const string_t& na)
{
  if(X1.strucType()!=_scalar || X2.strucType()!=_scalar)
    {
      where("SuTermVector::SuTermVector(Unknown, SuTermVector,SuTermVector)");
      error("term_nonscalar");
    }
  number_t nbc=vu.nbOfComponents();
  if(nbc<2)
    {
      where("SuTermVector::SuTermVector(Unknown, SuTermVector,SuTermVector)");
      error("unknown_vector",2);
    }
  if(X1.space_p!=X2.space_p)
    {
      where("SuTermVector::SuTermVector(Unknown, SuTermVector,SuTermVector)");
      error("term_mismatch_spaces", X1.space_p->name(), X2.space_p->name());
    }
  name_ = na;
  ValueType vt=X1.valueType();
  if(X2.valueType()==_complex) vt = _complex;
  space_p=X1.space_p;
  sulf_p=0;
  u_p=&vu;
  scalar_entries_p = 0;
  entries_p=new VectorEntry(vt, _vector, X1.nbDofs(), nbc);
  if(vt==_real)
    {
      Vector<real_t>::const_iterator it1 = X1.entries_p->rEntries_p->begin(),
                                     it2 = X2.entries_p->rEntries_p->begin();
      Vector<Vector<real_t> >::iterator itv = entries_p->rvEntries_p->begin();
      for(; itv!=entries_p->rvEntries_p->end(); ++itv,++it1,++it2)
        {(*itv)[0]=*it1; (*itv)[1]=*it2;}
    }
  else
    {
      Vector<Vector<complex_t> >::iterator itv = entries_p->cvEntries_p->begin();
      if(X1.valueType()==_real) // X2 is complex
        {
          Vector<real_t>::const_iterator it1 = X1.entries_p->rEntries_p->begin();
          Vector<complex_t>::const_iterator it2 = X2.entries_p->cEntries_p->begin();
          for(; itv!=entries_p->cvEntries_p->end(); ++itv,++it1,++it2)
            {(*itv)[0]=*it1; (*itv)[1]=*it2;}
        }
      else //X1 is complex
        {
          if(X2.valueType()==_real) // X2 is real
            {
              Vector<complex_t>::const_iterator it1 = X1.entries_p->cEntries_p->begin();
              Vector<real_t>::const_iterator it2 = X2.entries_p->rEntries_p->begin();
              for(; itv!=entries_p->cvEntries_p->end(); ++itv,++it1,++it2)
                {(*itv)[0]=*it1; (*itv)[1]=*it2;}
            }
          else //X3 complex
            {
              Vector<complex_t>::const_iterator it1 = X1.entries_p->cEntries_p->begin();
              Vector<complex_t>::const_iterator it2 = X2.entries_p->cEntries_p->begin();
              for(; itv!=entries_p->cvEntries_p->end(); ++itv,++it1,++it2)
                {(*itv)[0]=*it1; (*itv)[1]=*it2;}
            }
        }
    }
  computingInfo_.isComputed = true;
}

/*! construct vector SutermVector from n scalar SuTermVector's
    scalar SuTermVector's have to be defined on the same space
    vu has to be a vector unknown at least n */
SuTermVector::SuTermVector(const Unknown& vu, const std::list<const SuTermVector*>& Xs, const string_t& na)
{
  number_t nbX=Xs.size();
  if(nbX<2)
    {
      where("SuTermVector::SuTermVector(Unknown,list<const SuTermVector*>)");
      error("term_not_enough",2);
    }
  number_t nbc=vu.nbOfComponents();
  if(nbc<nbX)
    {
      where("SuTermVector::SuTermVector(Unknown,list<const SuTermVector*>)");
      error("unknown_vector",nbX);
    }
  Space* sp=0;
  ValueType vt=_real;
  std::list<const SuTermVector*>::const_iterator itX=Xs.begin();
  sp=(*itX)->space_p;
  number_t n=(*itX)->nbDofs();
  for(; itX!=Xs.end(); ++itX)
    {
      const SuTermVector& X=**itX;
      if(X.space_p!=sp)
        {
          where("SuTermVector::SuTermVector(Unknown,list<const SuTermVector*>)");
          error("term_mismatch_spaces", X.space_p->name(), sp->name());
        }
      if(X.valueType()==_complex) vt=_complex;
    }
  name_ = na;
  space_p=sp;
  sulf_p=0;
  u_p=&vu;
  scalar_entries_p = 0;
  entries_p=new VectorEntry(vt, _vector, n, nbc);
  if(vt==_real)  //all SuTermVector are real
    {
      std::list<Vector<real_t>::const_iterator> itlist;
      for(itX=Xs.begin(); itX!=Xs.end(); ++itX) itlist.push_back((*itX)->entries_p->rEntries_p->begin());
      std::list<Vector<real_t>::const_iterator>::iterator itl;
      Vector<Vector<real_t> >::iterator itv = entries_p->rvEntries_p->begin();
      for(; itv!=entries_p->rvEntries_p->end(); ++itv)
        {
          Vector<real_t>::iterator it=itv->begin();
          for(itl=itlist.begin(); itl!=itlist.end(); ++it,++itl) {*it=**itl; (*itl)++;}
        }
    }
  else // complex result
    {
      Vector<Vector<complex_t> >::iterator itv;
      number_t k=0;
      for(itX=Xs.begin(); itX!=Xs.end(); ++itX, k++)
        {
          const SuTermVector& X=**itX;
          if(X.valueType()==_real)
            {
              Vector<real_t>::const_iterator it=X.entries_p->rEntries_p->begin();
              for(itv = entries_p->cvEntries_p->begin(); itv!=entries_p->cvEntries_p->end(); ++itv,++it)(*itv)[k]=*it;
            }

          else
            {
              Vector<complex_t>::const_iterator it=X.entries_p->cEntries_p->begin();
              for(itv = entries_p->cvEntries_p->begin(); itv!=entries_p->cvEntries_p->end(); ++itv,++it)(*itv)[k]=*it;
            }
        }
    }
  computingInfo_.isComputed = true;
}

void SuTermVector::initFromFunction(const Unknown& u, const GeomDomain& dom, const OperatorOnFunction& opf,
                                    const string_t& na, bool noass)
{
  termType_ = _sutermVector;
  computingInfo_.noAssembly = noass;
  name_ = na;
  sulf_p = 0;
  u_p = &u;

  //find subspace or create it
  Space* spu = u.space();
  space_p = spu->findSubSpace(&dom, spu);
  if(space_p == 0) space_p = new Space(dom, *spu, spu->name() + "_" + dom.name());

  //allocate entries
  number_t n = space_p->nbDofs(), nv = number_t(u.nbOfComponents());
  entries_p = 0;
  scalar_entries_p = 0;
  if(n == 0) return;   // no entries

  //check structure of function
  if(opf.strucType()==_matrix)
    {
      where("SuTermVector::initFromFunction(Unknown, Domain, OperatorOnFunction, String, bool)");
      error("scalar_or_vector");
    }
  if(nv != opf.dims().first)
    {
      where("SuTermVector::initFromFunction(Unknown, Domain, OperatorOnFunction, String, bool)");
      error("operator_fct_mismatch_dims", nv, opf.dims().first);
    }

  //compute normal vectors if function required them (function parameters contains "_n" parameter)
  bool nor = opf.normalRequired(), norf = opf.funp()->requireNx;
  const MeshDomain* mdom = dom.meshDomain();
  if(mdom == 0)
    {
      where("SuTermVector::initFromFunction(Unknown, Domain, OperatorOnFunction, String, bool)");
      error("domain_notmesh", dom.name(), words("domain type", dom.domType()));
    }

  if(nor && !(mdom->isSideDomain() || mdom->spaceDim() > mdom->dim()))    //when function has "_n" as parameter but computation does not occur on a side domain
    {
      where("SuTermVector::initFromFunction(Unknown, Domain, OperatorOnFunction, String, bool)");
      warning("normal_ignored",mdom->name());
      nor = false;
    }
  std::vector<Vector<real_t> > ns;
  if(nor)
    {
      if(dom.meshDomain() == 0)
        {
          where("SuTermVector::initFromFunction(const Unknown&, const GeomDomain&, const OperatorOnFunction&, const String&, bool)");
          error("domain_notmesh", dom.name(), words("domain type", dom.domType()));
        }
      if(!dom.meshDomain()->orientationComputed) dom.meshDomain()->setNormalOrientation();
      interpolatedNormals(*space_p, ns);
    }

  //flag to transmit dof pointer tp function parameters
  bool assDof = opf.funp()->requireDof;

  //fill entries with f values
  Vector<real_t>* np = 0;
  ValueType vt = opf.valueType();
  StrucType st = opf.strucType();
  entries_p = new VectorEntry(vt, st, n, nv);
  if(vt == _real)
    {
      if(st == _scalar)
        {
          Vector<real_t>::iterator ite = entries_p->rEntries_p->begin();
#ifdef XLIFEPP_WITH_OMP
          #pragma omp parallel for firstprivate(np)
#endif // XLIFEPP_WITH_OMP
          for(number_t i = 0; i < n; i++)
            {
              if(nor)
                {
                  np = &ns[i];
                  if(norf) setNx(np);
                }
              if(assDof) setDof(const_cast<Dof*>(&space_p->dof(i + 1)));
              opf.eval(space_p->dof(i + 1).coords(), *(ite+i), np);
            }
        }
      if(st == _vector)
        {
          Vector<Vector<real_t> >::iterator ite = entries_p->rvEntries_p->begin();
#ifdef XLIFEPP_WITH_OMP
          #pragma omp parallel for firstprivate(np)
#endif // XLIFEPP_WITH_OMP
          for(number_t i = 0; i < n; i++)
            {
              if(nor)
                {
                  np = &ns[i];
                  if(norf) setNx(np);
                }
              if(assDof) setDof(const_cast<Dof*>(&space_p->dof(i + 1)));
              opf.eval(space_p->dof(i + 1).coords(), *(ite+i), np);
            }
        }
    }
  if(vt == _complex)
    {
      if(st == _scalar)
        {
          Vector<complex_t>::iterator ite = entries_p->cEntries_p->begin();
#ifdef XLIFEPP_WITH_OMP
          #pragma omp parallel for firstprivate(np)
#endif // XLIFEPP_WITH_OMP
          for(number_t i = 0; i < n; i++)
            {
              if(nor)
                {
                  np = &ns[i];
                  if(norf) setNx(np);
                }
              if(assDof) setDof(const_cast<Dof*>(&space_p->dof(i + 1)));
              opf.eval(space_p->dof(i + 1).coords(), *(ite+i), np);
            }
        }
      if(st == _vector)
        {
          Vector<Vector<complex_t> >::iterator ite = entries_p->cvEntries_p->begin();
#ifdef XLIFEPP_WITH_OMP
          #pragma omp parallel for firstprivate(np)
#endif // XLIFEPP_WITH_OMP
          for(number_t i = 0; i < n; i++)
            {
              if(nor)
                {
                  np = &ns[i];
                  if(norf) setNx(np);
                }
              if(assDof) setDof(const_cast<Dof*>(&space_p->dof(i + 1)));
              opf.eval(space_p->dof(i + 1).coords(), *(ite+i), np);
            }
        }
    }
  computingInfo_.isComputed = true;
}

/*! SuTermVector from the restriction of a SuTermVector to a given domain dom
    raise an error if dom does not intersect the SutermVector domain
    SutermVector sutv has to be computed and its entries_p has to be non zero
    In case of entries_p==0 and scalar_entries_p!=0 goto vector representation before calling this
*/
SuTermVector::SuTermVector(const SuTermVector& sutv, const GeomDomain& dom, const string_t& na)
{
  termType_ = _sutermVector;
  sulf_p = 0;
  u_p = 0;
  space_p = 0;
  entries_p = 0;
  scalar_entries_p = 0;
  if(!sutv.computed())
    {
      where("SuTermVector::SuTermVector");
      error("not_computed_term", sutv.name());
    }
  if(domain() == &dom)  //same domain
    {
      copy(sutv);
      where("SuTermVector::SuTermVector");
      warning("free_warning", "SuTermVector(SuTermVector,GeomDomain): no restriction to do, full copy");
    }

  Space* spu = sutv.u_p->space();
  space_p = spu->findSubSpace(&dom, spu);
  if(space_p == 0) space_p = new Space(dom, *spu, spu->name() + "_" + dom.name());
  u_p = sutv.u_p;
  number_t n = space_p->nbDofs(), nv = number_t(sutv.u_p->nbOfComponents());
  if(n == 0)
    {
      where("SuTermVector::SuTermVector");
      error("space_no_dof", space_p->name());
    }
  VectorEntry* ove = sutv.entries_p;
  VectorEntry* sove = sutv.scalar_entries_p;
  if(ove == 0)  //no standard entries
    {
      if(sove == 0)
        {
          where("SuTermVector::SuTermVector");
          error("term_no_entries");
        }
      else
        {
          where("SuTermVector::SuTermVector");
          error("entry_inconsistent_structure");
        }
    }

  std::vector<number_t> r = renumber(spu, space_p);
  number_t p = number_t(r.size());
  if(n == 0)  //no renumbering, copy all entries
    {
      entries_p = new VectorEntry(*sutv.entries_p);
    }
  else
    {
      entries_p = new VectorEntry(sutv.entries_p->valueType_, sutv.entries_p->strucType_, p, nv);
      std::vector<number_t>::iterator itr = r.begin();
      number_t k = 1;
      for(; itr != r.end(); itr++, k++)
        {
          number_t i = *itr;
          if(i != 0) entries_p->setValue(k, ove->getValue(i));
        }
    }
  computingInfo_.isComputed = true;
}

// SuTermVector from a linear combination of SuTermVector's
SuTermVector::SuTermVector(const LcTerm<SuTermVector>& lc)
{
  termType_ = _sutermVector;
  sulf_p = 0;
  u_p = 0;
  space_p = 0;
  entries_p = 0;
  scalar_entries_p = 0;
  compute(lc);
}

// assign operator (full copy)
SuTermVector& SuTermVector::operator=(const SuTermVector& suV)
{
  clear();
  copy(suV);
  return *this;
}

/*! construct SuTermVector from a SuTermVector with geometrical mapping, the unknowns may differ
    this process involves mesh interpolation and is restricted to Lagrange interpolation spaces
    it is assumed that the two unknowns have same structure
    - if fmap =0, identity is assumed
    - if errorOnOutDom = true the program stop on error, else the interpolated value is set to 0
*/
SuTermVector::SuTermVector(const Unknown& u, const GeomDomain& dom, const SuTermVector& sut,
                           const Function* fmap, bool errorOnOutDom, const string_t& na)
{
  termType_ = _sutermVector;
  name_ = na;
  sulf_p = 0;
  u_p = &u;
  Space* spu = u.space();
  space_p = spu->findSubSpace(&dom, spu);
  if(space_p == 0) space_p = new Space(dom, *spu, spu->name() + "_" + dom.name());
  number_t n = space_p->nbDofs(), nv = number_t(u.nbOfComponents());
  entries_p = 0;
  scalar_entries_p = 0;
  if(!sut.computed())
    {
      where("SuTermVector::SuTermVector");
      error("not_computed_term", "SuTermVector");
    }
  //if (nv!=sut.nbOfComponents())
  if(nv != sut.entries()->nbOfComponents)
    {
      where("SuTermVector::SuTermVector");
      error("term_mismatch_nb_components", nv, sut.entries()->nbOfComponents);
    }
  if(!spu->isFE())
    {
      where("SuTermVector::SuTermVector");
      error("not_fe_space_type", spu->name());
    }
  if(spu->feSpace()->interpolation()->type != _Lagrange)
    {
      where("SuTermVector::SuTermVector");
      error("lagrange_fe_space_only", spu->name());
    }
  ValueType vt = sut.valueType();
  StrucType st = sut.strucType();
  entries_p = new VectorEntry(vt, st, n, nv);
  //prepare interpolation
  Space* spdom = sut.space_p;
  const MeshDomain* mdom = spdom->domain()->meshDomain();
  mdom->buildKdTree();            //if not exist, create kdtree for fast location of points
  std::map<GeomElement*, const Element*> gelt2elt;

  for(number_t k = 0; k < spdom->nbOfElements(); k++)
    {
      const Element* elt = spdom->element_p(k);
      gelt2elt[elt->geomElt_p] = elt;
    }
  std::vector<real_t> v;
  //do interpolation
  if(vt == _real)
    {
      if(st == _scalar)
        {
          Vector<real_t>::iterator ite = entries_p->rEntries_p->begin();
          for(number_t i = 0; i < n; i++, ite++)
            {
              Point pt = space_p->dof(i + 1).coords();
              if(fmap != 0) pt = (*fmap)(space_p->dof(i + 1).coords(), v);  //map to dom
              GeomElement* gelt = mdom->locate(pt);
              if(gelt == 0)
                {
                  if(errorOnOutDom) locateError("in SutermVector constructor with mapping", *mdom, pt);
                  *ite = 0.;
                }
              else
                {
                  const Element* elt = gelt2elt[gelt];
                  elt->interpolate(*sut.entries_p, pt, elt->dofNumbers, *ite);
                }
            }
        }
      if(st == _vector)
        {
          Vector<Vector<real_t> >::iterator ite = entries_p->rvEntries_p->begin();
          Vector<real_t> zero(nv, 0.);
          for(number_t i = 0; i < n; i++, ite++)
            {
              Point pt = space_p->dof(i + 1).coords();
              if(fmap != 0) pt = (*fmap)(space_p->dof(i + 1).coords(), v);  //map to dom
              GeomElement* gelt = mdom->locate(pt);
              if(gelt == 0)
                {
                  if(errorOnOutDom) locateError("in SutermVector constructor with mapping", *mdom, pt);
                  *ite = zero;
                }
              else
                {
                  const Element* elt = gelt2elt[gelt];
                  elt->interpolate(*sut.entries_p, pt, elt->dofNumbers, *ite);
                }
            }
        }
    }
  if(vt == _complex)
    {
      if(st == _scalar)
        {
          Vector<complex_t>::iterator ite = entries_p->cEntries_p->begin();
          for(number_t i = 0; i < n; i++, ite++)
            {
              Point pt = space_p->dof(i + 1).coords();
              if(fmap != 0) pt = (*fmap)(space_p->dof(i + 1).coords(), v);  //map to dom
              GeomElement* gelt = mdom->locate(pt);
              if(gelt == 0)
                {
                  if(errorOnOutDom) locateError("in SutermVector constructor with mapping", *mdom, pt);
                  *ite = 0.;
                }
              else
                {
                  const Element* elt = gelt2elt[gelt];
                  elt->interpolate(*sut.entries_p, pt, elt->dofNumbers, *ite);
                }
            }
        }
      if(st == _vector)
        {
          Vector<Vector<complex_t> >::iterator ite = entries_p->cvEntries_p->begin();
          Vector<complex_t> zero(nv, 0.);
          for(number_t i = 0; i < n; i++, ite++)
            {
              Point pt = space_p->dof(i + 1).coords();
              if(fmap != 0) pt = (*fmap)(space_p->dof(i + 1).coords(), v);  //map to dom
              GeomElement* gelt = mdom->locate(pt);
              if(gelt == 0)
                {
                  if(errorOnOutDom) locateError("in SutermVector constructor with mapping", *mdom, pt);
                  *ite = zero;
                }
              else
                {
                  const Element* elt = gelt2elt[gelt];
                  elt->interpolate(*sut.entries_p, pt, elt->dofNumbers, *ite);
                }
            }
        }
    }
  computingInfo_.isComputed = true;
}

/*! produce the mapped SuTermVector from current GeomDomains to GeomDomain dom
    current SuTermVector has to be already computed
*/
SuTermVector* SuTermVector::mapTo(const GeomDomain& dom, const Unknown& u, bool errOnOut) const
{
  trace_p->push("SuTermVector::mapTo");
  if(!computed()) error("not_computed_term", "SuTermVector");
  SuTermVector* sut = 0;
  const GeomDomain* curdom = domain();
  if(curdom == 0) error("null_pointer", "domain");
  const Function* fmap = findMap(dom, *curdom);          // find map between domains
  //if (fmap==0) error("domain_map_undefined",dom.name(),curdom->name());
  //const GeomDomain* udom = u.space()->domain();
  //if (udom!= &dom && !udom->include(dom)) error("domain_not_subdomain",dom.name(),udom->name());
  //create new SuTermVector
  sut = new SuTermVector(u, dom, *this, fmap, errOnOut, name_ + "mapto_" + dom.name());
  trace_p->pop();
  return sut;
}

/*! extract unknown component term vector as a SuTermVector, create a  SuTermVector
    this function must be called from TermVector class in order to encapsulate SuTermVector in a TermVector
*/
SuTermVector SuTermVector::operator()(const ComponentOfUnknown& cu) const
{
  if(cu.parent() != u_p)
    {
      where("SuTermVector::operator()");
      error("term_wrong_unknown", name_, u_p->name(), cu.parent()->name());
    }
  if(cu.componentIndex() > nbOfComponents())
    {
      where("SuTermVector::operator()");
      error("term_wrong_unknown_component", name_, cu.componentIndex(), nbOfComponents());
    }
  if(nbOfComponents() == 1)
    {
      where("SuTermVector::operator()");
      error("unknown_not_scalar");
    }
  string_t na = name() + "_" + tostring(cu.componentIndex());
  SuTermVector su(na, &cu, space_p, valueType(), size(), 1, false); //create a scalar SuTermVector
  su.entries_p->extractComponent(*entries_p, cu.componentIndex());
  su.computed() = true;
  return su;
}

/*! move SuTermVector unknown to an other one
    consist in changing the unknown u of SuTermVector by the unknown v (u and v must be defined on the same Space!)
    with an optionnal management of components of unknown, examples
      u a scalar unknown and v a scalar unknown :  SuTermVector if not modified
      u a scalar unknown and v a 3 vector unknown, ci ={2} the SuTermVector is changed to a 3 column vector :
             new sutv=[ [0...0] [sutv] [0...0] ]
      u a 3-vector unknown and v a 3-vector unknown, ci ={2 3 1} the SuTermVector is changed columns are permuted :
             new sutv=[ [sutv2] [sutv3] [sutv1] ]   v may be equal to u in that case
      u a 4-vector unknown and v a 2-vector unknown, ci ={2 4} the SuTermVector is changed  :
             new sutv=[ [sutv2] [sutv4] ]
      u a 4-vector unknown and v a 3-vector unknown, ci ={2 0 4} the SuTermVector is changed :
             new sutv=[ [sutv2] [0...0] [sutv4] ]
    u is replaced by v,
    v may be equal to u but be cautious if the result is not of the dimension of v because many operations are not supported
    in doubt, create a new unknown of the consistent dimension
*/

void SuTermVector::changeUnknown(const Unknown& v, const Vector<number_t>& ci)
{
  dimen_t du = u_p->nbOfComponents(), dv = v.nbOfComponents();
  if(ci.size() == 0 || (du == 1 && dv == 1))
    {
      if(u_p != &v) u_p = &v;
      return;
    }
  if(entries_p == 0)
    {
      where("SuTermVector::changeUnknown");
      error("null_pointer", "entries_p");
    }
  entries_p->moveColumns(ci);
}

//destructor
SuTermVector::~SuTermVector()
{
  if(entries_p != 0) delete entries_p;       //delete entries
  if(scalar_entries_p != entries_p && scalar_entries_p != 0) delete scalar_entries_p;     //delete scalar entries
}

// size of SuTermvector counted in scalar
number_t SuTermVector::size() const
{
  if(entries_p != 0) return entries_p->size() * (u_p->nbOfComponents());
  if(scalar_entries_p != 0) return scalar_entries_p->size();
  return 0;
}


// size of Termvector counted in dofs
number_t SuTermVector::nbDofs() const
{
  return space_p->nbDofs();
}

//return n-th dof of unknown (n>=1)
const Dof& SuTermVector::dof(number_t n) const
{
  return space_p->dof(n);
}

// return geometric domain pointer attached to SuTermVector (may be 0)
const GeomDomain* SuTermVector::domain() const
{
  if(space_p == 0) return 0;  //no space -> no domain
  return space_p->domain();
}

// restriction of SuTermVector to a domain, return a pointer. May be null if dom does not intersect SuTermVector domain
SuTermVector* SuTermVector::onDomain(const GeomDomain& dom) const
{
  SuTermVector* sutv = new SuTermVector(*this, dom);
  if(sutv->size() == 0)
    {
      delete sutv;
      return 0;
    }
  else return sutv;
}

//deallocate memory used by a SuTermVector
void SuTermVector::clear()
{
  if(entries_p != 0) delete entries_p;
  if(scalar_entries_p != entries_p && scalar_entries_p != 0) delete scalar_entries_p;
  cdofs_.clear();
  entries_p = 0;
  scalar_entries_p = 0;
  computed() = false;
}

//! set of unknown spaces involved
std::set<const Space*> SuTermVector::unknownSpaces() const
{
  std::set<const Space*>sps;
  if(u_p != 0) sps.insert(u_p->space());
  return sps;
}

// return the actual pointer to entries (priority to scalar entry), return 0 is no allocated pointer
VectorEntry* SuTermVector::actual_entries() const
{
  VectorEntry* ve = scalar_entries_p;
  if(ve != 0) return ve;
  return entries_p;
}

//direct value managment (n>=1)
Value SuTermVector::getValue(number_t n) const
{
  if(entries_p == 0)
    {
      where("SuTermVector::getValue");
      error("null_pointer", "entries_p");
    }
  return entries_p->getValue(n);
}

void SuTermVector::setValue(number_t n, const Value& val)
{
  if(entries_p == 0)
    {
      where("SuTermVector::setValue");
      error("null_pointer", "entries_p");
    }
  entries_p->setValue(n, val);
}

/*! set dof values of SutermVector on a GeomDomainfrom from Value : u_i = v for all dof i on dom
    does not change other values
    dom has to be a subdomain of the domain of the currrent SuTermVector
    be cautious when Space is not a FE Lagrange space because does not imply that u=v on dom in the meaning of function*/
void SuTermVector::setValue(const Value& v, const GeomDomain& dom)
{
  //check space type
  Space* sp = u_p->space();

  //check value type
  StrucType st=strucType();
  ValueType vt=valueType();
  if(st!=v.strucType() || (vt==_real && v.valueType()!=_real))
    {
      where("SuTermVector::setValue(Value, Domain)");
      error("value_wrongtype", words("structure", st), "(" + words("value", vt) + ", " + words("value", v.valueType()) + ")",
            words("structure", v.strucType()), "(" + words("value", _real) + ", " + words("value", _real) + ")");
    }

  std::vector<number_t> rks;
  bool samedom= (domain() == &dom);
  if(!samedom)   //not same domain, get numberings
    {
      //get dof numbers of current SuTermVector, dof numbers of domain dom
      std::vector<number_t> curdofids = space_p->dofIds();
      std::vector<number_t> domdofids;
      Space* subsp=Space::findSubSpace(&dom,sp);
      if(subsp!=0) domdofids = subsp->dofIds();           //dofs of subspace related to dom
      else // built dof using dofsOn function
        {
          const_cast<GeomDomain&>(dom).updateParentOfSideElements();
          domdofids = sp->feSpace()->dofsOn(dom);        //dofs located on dom
        }
      number_t nd=domdofids.size();
      if(nd==0)
        {
          where("SuTermVector::setValue(Value, Domain)");
          error("dof_not_found");
        }
      rks.resize(nd,0);
      ranks(curdofids,domdofids,rks);   // costly find, to improve
    }

  //fill value
  VectorEntry* ve=entries_p;
  if(ve==0) ve=scalar_entries_p;
  if(ve==0)
    {
      where("SuTermVector::setValue(Value, Domain)");
      error("term_no_entries");
    }
  std::vector<number_t>::iterator itn=rks.begin();
  switch(st)
    {
      case _scalar:
        if(vt==_real)
          {
            real_t r = v.asReal();
            Vector<real_t>::iterator it= ve->rEntries_p->begin();
            if(samedom)
              for(; it!=ve->rEntries_p->end(); ++it) *it = r;
            else
              for(; itn!=rks.end(); ++itn)
                if(*itn!=0) *(it + *itn-1) = r;
          }
        else
          {
            complex_t c = 0.;
            if(v.valueType()==_real) c=v.asReal();
            else c=v.asComplex();
            Vector<complex_t>::iterator it= ve->cEntries_p->begin();
            if(samedom)
              for(; it!=ve->cEntries_p->end(); ++it) *it = c;
            else
              for(; itn!=rks.end(); ++itn)
                if(*itn!=0) *(it + *itn-1) = c;
          }
        break;
      case _vector:
        if(vt==_real)
          {
            Vector<real_t> rs=v.asRealVector();
            Vector<Vector<real_t> >::iterator it= ve->rvEntries_p->begin();
            if(samedom)
              for(; it!=ve->rvEntries_p->end(); ++it) *it = rs;
            else
              for(; itn!=rks.end(); ++itn)
                if(*itn!=0) *(it + *itn-1) = rs;
          }
        else
          {
            Vector<complex_t> cs;
            if(v.valueType()==_real) cs=v.asRealVector();
            else cs=v.asComplexVector();
            Vector<Vector<complex_t> >::iterator it= ve->cvEntries_p->begin();
            if(samedom)
              for(; it!=ve->cvEntries_p->end(); ++it) *it=cs;
            else
              for(; itn!=rks.end(); ++itn)
                if(*itn!=0) *(it + *itn-1) = cs;
          }
        break;
      default:
        where("SuTermVector::setValue(Value, Domain)");
        error("scalar_or_vector");
    }
}

/*! set dof values of GeomDomain from Function : u_i = f(M_i) for all dof i on dom, M_i point support of dof i
    does not change other values
    restricted to FE space
    dom has to be a subdomain of the domain of the currrent SuTermVector
    be cautious when Space is not a FE Lagrange space, because does not imply that u=f on dom in the meaning of function */
void SuTermVector::setValue(const Function& f, const GeomDomain& dom)
{
  //check space type
  Space* sp = u_p->space();
  if(sp->typeOfSpace() != _feSpace)
    {
      where("SuTermVector::setValue(Function, Domain)");
      error("not_fe_space_type", sp->name());
    }
  //check value type
  StrucType st=strucType();
  ValueType vt=valueType();
  if(st!=f.strucType() || (vt==_real && f.valueType()!=_real))
    {
      where("SuTermVector::setValue(Function, Domain)");
      error("value_wrongtype", words("structure", st), "(" + words("value", vt) + ", " + words("value", f.valueType()) + ")",
            words("structure", f.strucType()), "(" + words("value", _real) + ", " + words("value", _real) + ")");
    }

  std::vector<number_t> rks;
  bool samedom= (domain() == &dom);
  if(!samedom)   //not same domain, get numberings
    {
      //get dof numbers of current SuTermVector, dof numbers of domain dom
      std::vector<number_t> curdofids = space_p->dofIds();
      std::vector<number_t> domdofids;
      Space* subsp=Space::findSubSpace(&dom,sp);
      if(subsp!=0) domdofids = subsp->dofIds();           //dofs of subspace related to dom
      else // built dof using dofsOn function
        {
          const_cast<GeomDomain&>(dom).updateParentOfSideElements();
          domdofids = sp->feSpace()->dofsOn(dom);        //dofs located on dom
        }
      number_t nd=domdofids.size();
      if(nd==0)
        {
          where("SuTermVector::setValue(Function, Domain)");
          error("dof_not_found");
        }
      rks.resize(nd,0);
      ranks(curdofids,domdofids,rks);   // costly find, to improve
    }
  else rks=trivialNumbering<number_t>(1,space_p->nbDofs());

  //fill value
  VectorEntry* ve=entries_p;
  if(ve==0) ve=scalar_entries_p;
  if(ve==0)
    {
      where("SuTermVector::setValue(Function, Domain)");
      error("term_no_entries");
    }
  std::vector<number_t>::iterator itn=rks.begin();
  switch(st)
    {
      case _scalar:
        if(vt==_real)
          {
            Vector<real_t>::iterator it= ve->rEntries_p->begin();
            real_t r;
            for(; itn!=rks.end(); ++itn)
              if(*itn!=0) *(it + *itn-1) = f(sp->dof(*itn).coords(),r);
          }
        else
          {
            real_t r;
            complex_t c = 0.;
            Vector<complex_t>::iterator it= ve->cEntries_p->begin();
            for(; itn!=rks.end(); ++itn)
              if(*itn!=0)
                {
                  if(f.valueType()==_real) *(it + *itn-1) = f(sp->dof(*itn).coords(),r);
                  else *(it + *itn-1) = f(sp->dof(*itn).coords(),c);
                }
          }
        break;
      case _vector:
        if(vt==_real)
          {
            Vector<real_t> rs;
            Vector<Vector<real_t> >::iterator it= ve->rvEntries_p->begin();
            for(; itn!=rks.end(); ++itn)
              if(*itn!=0) *(it + *itn-1) = f(sp->dof(*itn).coords(),rs);
          }
        else
          {
            Vector<complex_t> cs;
            Vector<real_t> rs;
            Vector<Vector<complex_t> >::iterator it= ve->cvEntries_p->begin();
            for(; itn!=rks.end(); ++itn)
              if(*itn!=0)
                {
                  if(f.valueType()==_real) *(it + *itn-1) = f(sp->dof(*itn).coords(),rs);
                  else *(it + *itn-1) = f(sp->dof(*itn).coords(),cs);
                }
          }
        break;
      default:
        where("SuTermVector::setValue(Function, Domain)");
        error("scalar_or_vector");
    }
}

/*! set dof values of GeomDomain from SuTermVector : u_i = v_i for all dof i on dom
    does not change other values
    dom has to be a subdomain of the domain of the currrent SuTermVector
    be cautious when Space is not a FE Lagrange space, because does not imply that u=v on dom in the meaning of function*/
void SuTermVector::setValue(const SuTermVector& sut, const GeomDomain& dom)
{
  //check value type
  StrucType st=strucType();
  ValueType vt=valueType();
  if(st!=sut.strucType() || (vt==_real && sut.valueType()!=_real))
    {
      where("SuTermVector::setValue(SuTermVector, Domain)");
      error("value_wrongtype", words("structure", st), "(" + words("value", vt) + ", " + words("value", sut.valueType()) + ")",
            words("structure", sut.strucType()), "(" + words("value", _real) + ", " + words("value", _real) + ")");
    }

  //check space
  Space* sp = u_p->space();
  if(sut.u_p->space()!=sp)
    {
      where("SuTermVector::setValue(SuTermVector, Domain)");
      error("term_mismatch_spaces", sut.u_p->space()->name(), sp->name());
    }

  //build numberings
  std::vector<number_t> curdofids = space_p->dofIds();     //dof ids of current SuTermVector
  std::vector<number_t> sutdofids = sut.space_p->dofIds(); //dof ids of given SuTermVector
  std::vector<number_t> domdofids;                         //dof ids on domain
  Space* subsp=Space::findSubSpace(&dom,sp);
  if(subsp!=0) domdofids = subsp->dofIds();           //dofs of subspace related to dom
  else // built dof using dofsOn function
    {
      const_cast<GeomDomain&>(dom).updateParentOfSideElements();
      domdofids = sp->feSpace()->dofsOn(dom);        //dofs located on dom
    }
  std::vector<number_t> rks, rkc;
  ranks(sutdofids,domdofids,rks);   // ranks of domain dofs in sut dofs
  if(sut.domain()==&dom) rkc=trivialNumbering(number_t(1),sut.nbDofs());
  else ranks(curdofids,domdofids,rkc);   // ranks of domain dofs in current dofs

  //fill value
  VectorEntry* vec=entries_p;
  const VectorEntry* ves=sut.entries();
  if(vec==0) vec=scalar_entries_p;
  if(vec==0)
    {
      where("SuTermVector::setValue(SuTermVector, Domain)");
      error("term_no_entries");
    }
  if(ves==0) ves=scalar_entries_p;
  if(ves==0)
    {
      where("SuTermVector::setValue(SuTermVector, Domain)");
      error("term_no_entries");
    }
  std::vector<number_t>::iterator itns=rks.begin(), itnc=rkc.begin();
  switch(st)
    {
      case _scalar:
        if(vt==_real)
          {
            Vector<real_t>::iterator itc= vec->rEntries_p->begin();
            Vector<real_t>::const_iterator its= ves->rEntries_p->begin();
            for(; itns!=rks.end(); ++itns, ++itnc)
              if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
          }
        else
          {
            Vector<complex_t>::iterator itc= vec->cEntries_p->begin();
            if(sut.valueType()==_real)
              {
                Vector<real_t>::const_iterator its= ves->rEntries_p->begin();
                for(; itns!=rks.end(); ++itns, ++itnc)
                  if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
              }
            else
              {
                Vector<complex_t>::const_iterator its= ves->cEntries_p->begin();
                for(; itns!=rks.end(); ++itns, ++itnc)
                  if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
              }
          }
        break;
      case _vector:
        if(vt==_real)
          {
            Vector<Vector<real_t> >::iterator itc= vec->rvEntries_p->begin();
            Vector<Vector<real_t> >::const_iterator its= ves->rvEntries_p->begin();
            for(; itns!=rks.end(); ++itns, ++itnc)
              if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
          }
        else
          {
            Vector<Vector<complex_t> >::iterator itc= vec->cvEntries_p->begin();
            if(sut.valueType()==_real)
              {
                Vector<Vector<real_t> >::const_iterator its= ves->rvEntries_p->begin();
                for(; itns!=rks.end(); ++itns, ++itnc)
                  if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
              }
            else
              {
                Vector<Vector<complex_t> >::const_iterator its= ves->cvEntries_p->begin();
                for(; itns!=rks.end(); ++itns, ++itnc)
                  if(*itns!=0 && *itnc!=0) *(itc + *itnc-1) = *(its+*itns-1);
              }
          }
        break;
      default:
        where("SuTermVector::setValue(SuTermVector, Domain)");
        error("scalar_or_vector");
    }
}

/*! set dof values of SutermVector from a given SuTermVector : u_i = v_i for all dof i of the given SutermVector
    does not change other values */
void SuTermVector::setValue(const SuTermVector& sut)
{
  //check value type
  StrucType st=strucType();
  ValueType vt=valueType();
  if(st!=sut.strucType() || (vt==_real && sut.valueType()!=_real))
    {
      where("SuTermVector::setValue(SuTermVector)");
      error("value_wrongtype", words("structure", st), "(" + words("value", vt) + ", " + words("value", sut.valueType()) + ")",
            words("structure", sut.strucType()), "(" + words("value", _real) + ", " + words("value", _real) + ")");
    }

  //check space
  Space* sp = u_p->space();
  if(sut.u_p->space()!=sp)
    {
      where("SuTermVector::setValue(SuTermVector)");
      error("term_mismatch_spaces", sut.u_p->space()->name(), sp->name());
    }

  //build numberings
  std::vector<number_t> rks;
  if(sp==space_p) rks = sut.space_p->dofIds();  //dof ids of given SuTermVector
  else
    {
      std::vector<number_t> curdofids = space_p->dofIds();     //dof ids of current SuTermVector
      std::vector<number_t> sutdofids = sut.space_p->dofIds(); //dof ids of given SuTermVector
      ranks(curdofids,sutdofids,rks);   // ranks of sut dofs in current dofs
    }

  //fill value
  VectorEntry* vec=entries_p;
  const VectorEntry* ves=sut.entries();
  if(vec==0) vec=scalar_entries_p;
  if(vec==0)
    {
      where("SuTermVector::setValue(SuTermVector)");
      error("term_no_entries");
    }
  if(ves==0) ves=scalar_entries_p;
  if(ves==0)
    {
      where("SuTermVector::setValue(SuTermVector)");
      error("term_no_entries");
    }
  std::vector<number_t>::iterator itn=rks.begin();
  switch(st)
    {
      case _scalar:
        if(vt==_real)
          {
            Vector<real_t>::iterator itc= vec->rEntries_p->begin();
            Vector<real_t>::const_iterator its= ves->rEntries_p->begin();
            for(; itn!=rks.end(); ++itn, ++its)
              if(*itn!=0) *(itc + *itn-1) = *its;
          }
        else
          {
            Vector<complex_t>::iterator itc= vec->cEntries_p->begin();
            if(sut.valueType()==_real)
              {
                Vector<real_t>::const_iterator its= ves->rEntries_p->begin();
                for(; itn!=rks.end(); ++itn, ++its)
                  if(*itn!=0) *(itc + *itn-1) = *its;
              }
            else
              {
                Vector<complex_t>::const_iterator its= ves->cEntries_p->begin();
                for(; itn!=rks.end(); ++itn, ++its)
                  if(*itn!=0) *(itc + *itn-1) = *its;
              }
          }
        break;
      case _vector:
        if(vt==_real)
          {
            Vector<Vector<real_t> >::iterator itc= vec->rvEntries_p->begin();
            Vector<Vector<real_t> >::const_iterator its= ves->rvEntries_p->begin();
            for(; itn!=rks.end(); ++itn, ++its)
              if(*itn!=0) *(itc + *itn-1) = *its;
          }
        else
          {
            Vector<Vector<complex_t> >::iterator itc= vec->cvEntries_p->begin();
            if(sut.valueType()==_real)
              {
                Vector<Vector<real_t> >::const_iterator its= ves->rvEntries_p->begin();
                for(; itn!=rks.end(); ++itn, ++its)
                  if(*itn!=0) *(itc + *itn-1) = *its;
              }
            else
              {
                Vector<Vector<complex_t> >::const_iterator its= ves->cvEntries_p->begin();
                for(; itn!=rks.end(); ++itn, ++its)
                  if(*itn!=0) *(itc + *itn-1) = *its;
              }
          }
        break;
      default:
        where("SuTermVector::setValue(SuTermVector)");
        error("scalar_or_vector");
    }
}


// evaluate interpolated value at a point, use expansive interpolation method
Value SuTermVector::evaluate(const Point& p) const
{
  if(entries_p == 0)
    {
      where("SuTermVector::evaluate");
      error("null_pointer", "entries_p");
    }
  const FeSpace* fsp = space_p->feSpace();
  const FeSubSpace* ssp = space_p->feSubSpace();
  if(fsp == 0 && ssp == 0)
    {
      where("SuTermVector::evaluate");
      error("not_fe_space_type", space_p->name());
    }
  dimen_t dimf, dimu;
  if(fsp != 0) dimf = fsp->dimFun();
  else dimf = ssp->dimFun();
  if(dimf >1) dimf += space_p->dimPoint()-space_p->dimDomain();  //modify dimension because of mapping mD->nD
  dimu = u_p->nbOfComponents();
  dimen_t d = dimf * dimu;

  if(entries_p->rEntries_p != 0)  //real scalar input
    {
      if(d == 1)  //real scalar output
        {
          real_t r;
          if(fsp != 0) fsp->interpolate(*entries_p->rEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->rEntries_p, p, r, _id);
          return Value(r);
        }
      else  //real vector output
        {
          Vector<real_t> r(d);
          if(fsp != 0) fsp->interpolate(*entries_p->rEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->rEntries_p, p, r, _id);
          return Value(r);
        }
    }

  if(entries_p->cEntries_p != 0)  //complex scalar input
    {
      if(d == 1)  //complex scalar output
        {
          complex_t r;
          if(fsp != 0) fsp->interpolate(*entries_p->cEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->cEntries_p, p, r, _id);
          return Value(r);
        }
      else //complex vector output
        {
          Vector<complex_t> r(d);
          if(fsp != 0) fsp->interpolate(*entries_p->cEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->cEntries_p, p, r, _id);
          return Value(r);
        }
    }

  if(entries_p->rvEntries_p != 0)  //real vector input
    {
      if(d == 1)  //real scalar output
        {
          Vector<real_t> r(1);
          if(fsp != 0) fsp->interpolate(*entries_p->rvEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->rvEntries_p, p, r, _id);
          return Value(r[0]);
        }
      else //real vector output
        {
          Vector<real_t> r(d);
          if(fsp != 0) fsp->interpolate(*entries_p->rvEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->rvEntries_p, p, r, _id);
          return Value(Vector<real_t>(r));
        }
    }

  if(entries_p->cvEntries_p != 0)  //complex vector input
    {
      if(d == 1)  //complex scalar output
        {
          Vector<complex_t> r(1);
          if(fsp != 0) fsp->interpolate(*entries_p->cvEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->cvEntries_p, p, r, _id);
          return Value(r[0]);
        }
      else //complex vector output
        {
          Vector<complex_t> r(d);
          if(fsp != 0) fsp->interpolate(*entries_p->cvEntries_p, p, r, _id);
          else ssp->interpolate(*entries_p->cvEntries_p, p, r, _id);
          return Value(Vector<complex_t>(r));
        }
    }

  where("SuTermVector::evaluate");
  error("term_no_entries");
  return Value(0.);
}


/*!build the required subspaces and the largest subspace (possibly the whole space)
  this function :
   - travel along basic linear forms of the linear combination
   - identify the largest geometric domain
   - build the largest subspace (for SuTermVector dofs numbering)
   - create subspace related to each basic linear form if it does not exist
   - create dofs numbering between domain subspaces ands largest subspace

 dofs numbering :
    whole space     largest subspace    domain subspace
            1               5                2
            2               3                1
            3               1
            4
            5
*/
void SuTermVector::buildSubspaces()
{
  //for safety
  if(sulf_p == 0) return ;               //no linear form, nothing to do
  trace_p->push("SuTermVector::buildSubspaces");

  //case of not a FeSpace
  Space* spu = sulf_p->unknown()->space(); //space of the unknown
  number_t nblf = sulf_p->size();            //number of basic linear forms in linear form
  std::vector<const GeomDomain*> doms(nblf, 0);
  space_p = spu;
  if(spu->typeOfSpace() != _feSpace)
    {
      spu->buildSubspaces(doms, subspaces); //build subspaces
      trace_p->pop();
      return;
    }

  //only for FeSpace, make list of domains of each basic linear form
  number_t n = 0;
  for(it_vlfp it = sulf_p->begin(); it != sulf_p->end(); it++, n++)  //travel the linear combination of forms
    {
      BasicLinearForm* blc = it->first;
      switch(blc->type())
        {
          case _intg :
            {
              if(blc->asIntgForm()->computationType()==_FEextComputation && blc->asIntgForm()->opu()->extensionRequired())
                doms[n] = &blc->asIntgForm()->domain()->meshDomain()->extendDomain();  //find or create extension of domain dom
              else doms[n] = blc->asIntgForm()->domain();
              break;
            }
          //case _bilinearAsLinear :
          //doms[n] = blc->asBilinearAsLinearForm()->domain(); //get the domain of linear form
          // break;
          case _doubleIntg :                                   //to do later
          default          :
            error("lform_not_handled", words("form type", blc->type()));      //unknown basic linear form
        }
    }

  space_p = spu->buildSubspaces(doms, subspaces); //build subspaces

  trace_p->pop();
}

/*! if vector unknown, create vector representation from scalar representation (have to exist)
    if scalar unknown, set entries_p = scalar_entries_p
    keepScalar : if false (default value), scalar entries are deleted
*/
void SuTermVector::toVector(bool keepScalar)
{
  if(scalar_entries_p == 0) return;  // no scalar representation, nothing done
  trace_p->push("SuTermVector::toVector()");
  if(entries_p != 0 && scalar_entries_p != entries_p) delete entries_p;  // delete old entries
  dimen_t nbc = u_p->nbOfComponents();
  if(nbc <= 1)  //not a vector unknown
    {
      entries_p = scalar_entries_p;
      trace_p->pop();
      return;
    }
  //reallocate entries
  ValueType vt = scalar_entries_p->valueType_;
  std::vector<number_t> dofs = space_p->dofIds();
  entries_p = new VectorEntry(vt, nbc, dofs.size());
  std::map<number_t, number_t> dofindex;
  std::vector<number_t>::iterator itd = dofs.begin();
  number_t n = 1;
  for(; itd != dofs.end(); itd++, n++) dofindex[*itd] = n;
  std::vector<DofComponent>::iterator itcd = cdofs_.begin();
  std::map<number_t, number_t>::iterator itm;
  if(vt == _real)  //real case
    {
      Vector<real_t>& sval = *scalar_entries_p->rEntries_p;
      Vector<Vector<real_t> >& val = *entries_p->rvEntries_p;
      Vector<real_t>::iterator itv = sval.begin();
      for(; itcd != cdofs_.end(); itcd++, itv++)
        {
          itm = dofindex.find(itcd->dofnum);
          if(itm == dofindex.end())
            error("dof_not_found");
          val(itm->second)(itcd->numc) = *itv;
        }
    }
  else //complex case
    {
      Vector<complex_t>& sval = *scalar_entries_p->cEntries_p;
      Vector<Vector<complex_t> >& val = *entries_p->cvEntries_p;
      Vector<complex_t>::iterator itv = sval.begin();
      for(; itcd != cdofs_.end(); itcd++, itv++)
        {
          itm = dofindex.find(itcd->dofnum);
          if(itm == dofindex.end())
            error("dof_not_found");
          val(itm->second)(itcd->numc) = *itv;
        }
    }
  if(!keepScalar)
    {
      delete scalar_entries_p;
      scalar_entries_p = 0;
    }
  trace_p->pop();
}

/*! create scalar representation scalar_entries_p and cdofs_ numbering vectors
    - if unknown is a vector unknown a new VectorEntry is created using  VectorEntry::toScalar() function
      and linked to the  VectorEntry pointer scalar_entries_p. In that case, if keepVector is false the vector representation is destroyed
    - if unknowns is a scalar unknown, no new VectorEntry : scalar_entries_p = entries_p
    This function builds in any case, the cdofs_ numbering vector (vector of ComponentDof's)
*/
void SuTermVector::toScalar(bool keepVector)
{
  if(scalar_entries_p != 0) return;  //already in scalar form

  trace_p->push("SuTermVector::toScalar");
  //create or update scalar_entries_p
  dimen_t nbc = u_p->nbOfComponents();
  if(nbc == 1)  scalar_entries_p = entries_p;                       //no creation, pointer copy
  else if(entries_p != 0) scalar_entries_p = entries_p->toScalar();  //create new scalar entries from vector entries

  //create cdofs_u
  std::vector<number_t> dofs = space_p->dofIds();
  cdofs_.resize(dofs.size()*nbc);
  std::vector<number_t>::iterator itd = dofs.begin();
  std::vector<DofComponent>::iterator itdc = cdofs_.begin();
  for(; itd != dofs.end(); itd++)
    for(dimen_t d = 0; d < nbc; d++, itdc++) *itdc = DofComponent(u_p, *itd, d + 1);

  if(!keepVector && entries_p != scalar_entries_p)
    {
      delete entries_p;
      entries_p = 0;
    }
  trace_p->pop();
}

/*! when current SutermVector has a component unknown, return vector unknown representation
    if is not a component unknown SuTermVector, a copy of current SuTermVector is returned
*/
SuTermVector SuTermVector::toVectorUnknown() const
{
  trace_p->push("SuTermVector::toVectorUnknown");
  SuTermVector res(*this);
  if(!u_p->isComponent() || entries_p == 0)
    {
      trace_p->pop();
      return res;
    }
  res.entries_p->toVector(u_p->parent()->nbOfComponents(), u_p->componentIndex());
  res.u_p = u_p->parent();
  trace_p->pop();
  return res;
}

// save SuTermVector to file
void SuTermVector::saveToFile(const string_t& fn, bool encode) const
{
  if(entries_p != 0) entries_p->saveToFile(fn, encode);
  if(scalar_entries_p != 0) scalar_entries_p->saveToFile(fn, encode);
}

//return value type (_real,_complex)
ValueType SuTermVector::valueType() const
{
  if(entries_p != 0) return entries_p->valueType_;
  if(scalar_entries_p != 0) return scalar_entries_p->valueType_;
  if(sulf_p != 0) return sulf_p->valueType();
  return _none;
}
//return structure type (_scalar,_vector)
StrucType SuTermVector::strucType() const
{
  if(entries_p != 0) return entries_p->strucType_;
  if(scalar_entries_p != 0) return scalar_entries_p->strucType_;
  if(sulf_p != 0) return sulf_p->strucType();
  return _scalar;
}

/* ----------------------------------------------------------------------------------------
  compute abs, real and imaginary part of SuTermVector
  note : current SuTermVector is modified
  ----------------------------------------------------------------------------------------*/
SuTermVector& SuTermVector::toAbs()
{
  if(entries_p != 0) entries_p->toAbs();
  if(scalar_entries_p != 0) scalar_entries_p->toAbs();
  return *this;
}

SuTermVector& SuTermVector::toReal()
{
  if(entries_p != 0) entries_p->toReal();
  if(scalar_entries_p != 0) scalar_entries_p->toReal();
  return *this;
}

SuTermVector& SuTermVector::toImag()
{
  if(entries_p != 0) entries_p->toImag();
  if(scalar_entries_p != 0) scalar_entries_p->toImag();
  return *this;
}

SuTermVector& SuTermVector::toComplex()
{
  if(entries_p != 0) entries_p->toComplex();
  if(scalar_entries_p != 0) scalar_entries_p->toComplex();
  return *this;
}

// take the conjugate of entries, if no entries nothing is done
// when conjugate is taken, the linear form is not updated to take it into account !
SuTermVector& SuTermVector::toConj()
{
  if(valueType()==_real) return *this;  //nothing to do
  if(entries_p != 0) entries_p->toConj();
  if(scalar_entries_p != 0) scalar_entries_p->toConj();
  return *this;

}
// round to zero all coefficients close to 0 (|.| < aszero)
SuTermVector& SuTermVector::roundToZero(real_t aszero)
{
  if(entries_p != 0) entries_p->roundToZero(aszero);
  if(scalar_entries_p != 0) scalar_entries_p->roundToZero(aszero);
  return *this;

}

/* ----------------------------------------------------------------------------------------
  compute SuTermVector
  remember that it deal with a single unknown
  this routine first analyse the linear form to be computed,
  it collects computations by domain and type of computation

  !!! For the moment, only FE computation is handled (not the BEM and spectral methods)
  ----------------------------------------------------------------------------------------*/
void SuTermVector::compute()
{
  if(sulf_p == 0) return;     //no linear form to compute
  if(sulf_p->size() == 0) return;    //void SuTermVector !!!
  trace_p->push("SuTermVector::compute");
  //retry unknown
  const Unknown* up = sulf_p->unknown();

  //analyse the combination of linear forms to be computed and split them by domain (and by type of computation)
  std::map<const GeomDomain*, std::vector<lfPair> > intgByDomain;  //list of linear form by domain
  std::map<const GeomDomain*, std::vector<lfPair> >::iterator itm;
  std::list<lfPair>  bilinearAslinear;  //list of bilinear_as_linear forms
  std::list<lfPair>::iterator itl;
  ValueType vt = _real;
  StrucType st = sulf_p->begin()->first->strucType();
  dimPair dims = dimPair(0,0);
  for(it_vlfp it = sulf_p->begin(); it != sulf_p->end(); it++)
    {
      switch(it->first->type())
        {
          case _intg :  //integral type
            {
              IntgLinearForm* ifl = it->first->asIntgForm(); //downcast
              const GeomDomain* dom = ifl->domain();
              itm = intgByDomain.find(dom);
              if(itm == intgByDomain.end())
                intgByDomain[dom] = std::vector<lfPair>(1, lfPair(*it));
              else
                intgByDomain[dom].push_back(lfPair(*it));
              //value and structure types
              if(it->first->valueType() == _complex || it->second.imag() != 0) vt = _complex;
              if(dims.first==0) dims=ifl->opu()->dimsRes();
              else if(dims!=ifl->opu()->dimsRes())
                error("operator_mismatch_size", dims.first, ifl->opu()->dimsRes().first);
              if(it->first->strucType() != st)
                error("term_mismatch_structures", words("structure", it->first->strucType()), words("structure", st));
              break;
            }
          case _bilinearAsLinear :
            {
              bilinearAslinear.push_back(lfPair(*it));
              if(it->first->valueType() == _complex) vt = _complex;
              if(it->first->strucType() != st)
                error("term_mismatch_structures", words("structure", it->first->strucType()), words("structure", st));
              break;
            }
          default :
            error("lform_not_handled", words("form type", it->first->type()));

        }
    }
  std::vector<SuLinearForm> sulfs;   //list of SuLinearForm, each of them is defined on a unique domain
  for(itm = intgByDomain.begin(); itm != intgByDomain.end(); itm++)
    {
      sulfs.push_back(SuLinearForm(itm->second));
    }
  intgByDomain.clear();

  //create entries
  if(dims.second>1) error("scalar_or_vector");
  if(up->nbOfComponents()>1 && dims.first>1)
    error("term_inconsistent_nb_components");
  dimen_t du = std::max(up->nbOfComponents(),dims.first);
  if(du == 1)
    {
      st=_scalar;
      entries_p = new VectorEntry(sulf_p->valueType(), _scalar, space_p->dimSpace());
    }
  else
    {
      st=_vector;
      entries_p = new VectorEntry(sulf_p->valueType(), _vector, space_p->nbDofs(), du);
    }

  //computation of SuLinearForms from intg
  std::vector<SuLinearForm>::iterator itsulf;
  UnknownType ut = up->type();
  //st = up->strucType();
  real_t r = 0;
  complex_t c = 0.;
  for(itsulf = sulfs.begin(); itsulf != sulfs.end(); itsulf++)
    {
      switch(ut)
        {
          case _feUnknown :
            switch(st)
              {
                case _scalar :
                  if(vt == _real) computeFE(*itsulf, *(entries_p->rEntries_p), r);
                  else computeFE(*itsulf, *(entries_p->cEntries_p), c);
                  break;
                case _vector :
                  if(vt == _real) computeFE(*itsulf, *(entries_p->rvEntries_p), r);
                  else computeFE(*itsulf, *(entries_p->cvEntries_p), c);
                  break;
                default :
                  error("scalar_or_vector");
              }
            break;
          case _spUnknown :
          default :
            error("unknown_not_handled", words("unknown", ut));
        }
    }

  // //computation of SuLinearForms from bilinear_as_linear form
  // for (itl=bilinearAslinear.begin(); itl!=bilinearAslinear.end(); ++itl)
  //   switch(st)
  //     {
  //       case _scalar :
  //         if (vt == _real) computeBilAsLin(*itl, *(entries_p->rEntries_p), r);
  //         else computeBilAsLin(*itl, *(entries_p->cEntries_p), c);
  //         break;
  //       case _vector :
  //         if (vt == _real) computeBilAsLin(*itl, *(entries_p->rvEntries_p), r);
  //         else computeBilAsLin(*itl, *(entries_p->cvEntries_p), c);
  //         break;
  //       default :
  //         error("scalar_or_vector");
  //     }

  computed() = true;
  trace_p->pop();
}

/*! extend current vector to numbering of vector v (useful for algebraic computation on SuTermVector)
    update entries and not scalar entries if allocated
*/
void SuTermVector::extendTo(const SuTermVector& v)
{
  if(&v == this) return;             //same vector nothing to extend
  if(space_p == v.space_p) return;   //same space nothing to extend
  if(u_p != v.u_p && u_p != v.u_p->dual_p())
    {
      where("SuTermVector::extendTo");
      error("term_inconsistent_unknowns");
    }
  if(entries_p == 0) return;  //entries not allocated, do nothing

  //extension may be required
  std::vector<number_t> dofrenumber = renumber(v.space_p, space_p);
  if(dofrenumber.size() == 0) return;  //a trivial renumbering
  entries_p->extendEntries(dofrenumber, v.space_p->dimSpace());
  space_p = v.space_p;   //update SuTermVector space
}

/*! extend current vector to numbering of space sp
   update entries and not scalar entries
*/
void SuTermVector::extendTo(const Space& sp)
{
  if(space_p == &sp) return;   //same space nothing to extend
  if(space_p->rootSpace() != sp.rootSpace())
    error("term_mismatch_spaces", space_p->rootSpace()->name(), sp.rootSpace()->name());
  if(entries_p == 0) return;  //entries not allocated, do nothing

  std::vector<number_t> dofrenumber = renumber(&sp, space_p);
  if(dofrenumber.size() == 0) return;  //trivial renumbering
  entries_p->extendEntries(dofrenumber, sp.dimSpace());
  space_p = const_cast<Space*>(&sp);    //update SuTermVector space
}

/*! extend current vector to numbering given by DofComponent vector ndofs
    if useDual is true, use dual components dof (default = false)
    update scalar entries and not vector entries
    if scalar_entries = 0  move to scalar representation
    Note : differs from previous extendTo functions by the fact it works with component dof numbering instead of dof numbering
*/
void SuTermVector::extendScalarTo(const std::vector<DofComponent>& ndofs, bool useDual)
{
  VectorEntry* b = scalar_entries();
  if(b == 0)
    {
      toScalar();
      b = scalar_entries();
    }
  std::vector<number_t> renum;
  if(!useDual) renum = renumber(ndofs, cdofs());
  else         renum = renumber(ndofs, dualDofComponents(cdofs()));
  if(renum.size() != 0)
    {
      b->extendEntries(renum, ndofs.size());      //extend vector to matrix column numbering
      cdofs() = ndofs;       //update cdofs
    }
}

/*! adjust scalar entries to cdofs numbering given
    same as extendScalarTo but different method, see general adjustScalarEntries function
*/
void SuTermVector::adjustScalarEntries(const std::vector<DofComponent>& ncdofs)
{
  adjustScalarEntriesG(scalar_entries_p, cdofs_, ncdofs);
}

//algebraic operations on SuTermVectors
SuTermVector& SuTermVector::operator+=(const SuTermVector& v)           // operation U+=V
{
  if(u_p != v.u_p && u_p != v.u_p->dual_p())
    {
      where("SutermVector::operator+=");
      error("term_inconsistent_unknowns"); //different unknowns, error
    }
  if(v.entries_p == 0 && v.scalar_entries_p == 0) return *this;    //v has no entries, add 0 !

  // first, we have to check if both vectors have the same entries structure
  // if no, we have to set v to the same structure, through a copy
  bool copyNeeded=false;
  bool vToScalar=false;
  bool vToVector=false;
  if(scalar_entries_p != 0 && scalar_entries_p != entries_p && v.scalar_entries_p == 0)
    {
      copyNeeded=vToScalar=true;
    }
  else if(scalar_entries_p == 0 && v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
    {
      copyNeeded=vToVector=true;
    }

  if(copyNeeded)
    {
      SuTermVector vcopy(v);
      if(vToScalar)
        {
          vcopy.toScalar();
        }
      if(vToVector)
        {
          vcopy.toVector();
        }

      Space* spu = space_p, *spv = v.space_p;
      if(spu == spv)  // same numbering => go to entries operation
        {
          if(vToVector)
            {
              if(entries_p == 0) entries_p = new VectorEntry(vcopy.valueType(), vcopy.strucType(), vcopy.nbDofs(), vcopy.nbOfComponents());
              *entries_p += *vcopy.entries_p;
            }
          if(vToScalar)
            {
              if(scalar_entries_p == 0) scalar_entries_p = new VectorEntry(vcopy.valueType(), vcopy.strucType(), vcopy.nbDofs()*vcopy.nbOfComponents(), 1);
              *scalar_entries_p += *vcopy.scalar_entries_p;
            }
        }
      else // not the same numbering => dynamic renumbering
        {
          Space* spuv = mergeSubspaces(spu,spv); // create union of subspaces, may be spu or spv or anything else
          space_p=spuv;                          // update subspace member

          if(spu != spuv)  // extend current entry if required
            {
              if(vToVector)
                {
                  if(entries_p != 0 && vcopy.entries_p != 0)
                    {
                      std::vector<number_t> dofrenumber = renumber(spuv, spu);
                      if(dofrenumber.size() != 0) entries_p->extendEntries(dofrenumber, spuv->dimSpace());
                    }
                  VectorEntry* ve = vcopy.entries_p;

                  if(spv != spuv)
                    {
                      if(entries_p != 0 && vcopy.entries_p != 0)
                        {
                          ve = new VectorEntry(*vcopy.entries_p); // copy of v entries
                          std::vector<number_t> dofrenumber = renumber(spuv, spv);
                          if(dofrenumber.size() != 0) ve->extendEntries(dofrenumber, spuv->dimSpace());  // do algebraic operation on entries
                        }
                    }
                  if(entries_p != 0 && vcopy.entries_p != 0)
                    {
                      *entries_p += *ve;  // do algebraic operation on entries
                    }
                }
              if(vToScalar)
                {
                  if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                    {
                      std::vector<number_t> dofrenumber;
                      if(u_p == vcopy.u_p) dofrenumber=renumber(cdofs(), vcopy.cdofs());
                      else                  dofrenumber=renumber(cdofs(), dualDofComponents(vcopy.cdofs()));
                      if(dofrenumber.size() != 0) scalar_entries_p->extendEntries(dofrenumber, cdofs().size());
                    }
                  VectorEntry* vse = vcopy.scalar_entries_p;

                  if(spv != spuv)
                    {
                      if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                        {
                          vse=new VectorEntry(*vcopy.scalar_entries_p);
                          std::vector<number_t> dofrenumber;
                          if(u_p == vcopy.u_p) dofrenumber=renumber(vcopy.cdofs(), cdofs());
                          else                  dofrenumber=renumber(vcopy.cdofs(), dualDofComponents(cdofs()));
                          if(dofrenumber.size() != 0) vse->extendEntries(dofrenumber, vcopy.cdofs().size());
                        }
                    }
                  if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                    {
                      *scalar_entries_p += *vse;  // do algebraic operation on entries
                    }
                }
            }
        }
    }
  else
    {
      Space* spu = space_p, *spv = v.space_p;
      if(spu == spv)  // same numbering => go to entries operation
        {
          if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
            {
              if(scalar_entries_p == 0) scalar_entries_p = new VectorEntry(v.valueType(), v.strucType(), v.nbDofs()*v.nbOfComponents(), 1);
              if(scalar_entries_p != entries_p) *scalar_entries_p += *v.scalar_entries_p;
            }
          if(v.entries_p != 0)
            {
              if(entries_p == 0) entries_p = new VectorEntry(v.valueType(), v.strucType(), v.nbDofs(), v.nbOfComponents());
              *entries_p += *v.entries_p;
            }
        }
      else // not the same numbering => dynamic renumbering
        {
          Space* spuv = mergeSubspaces(spu, spv); // create union of subspaces, may be spu or spv or anything else
          space_p = spuv;                         // update subspace member

          if(spu != spuv)  //extend current entry if required
            {
              if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
                {
                  std::vector<number_t> dofrenumber;
                  if(u_p == v.u_p) dofrenumber=renumber(cdofs(), v.cdofs());
                  else              dofrenumber=renumber(cdofs(), dualDofComponents(v.cdofs()));
                  if(dofrenumber.size() != 0) scalar_entries_p->extendEntries(dofrenumber, cdofs().size());
                }
              if(v.entries_p != 0)
                {
                  std::vector<number_t> dofrenumber = renumber(spuv, spu);
                  if(dofrenumber.size() != 0) entries_p->extendEntries(dofrenumber, spuv->dimSpace());
                }
            }
          VectorEntry* ve = v.entries_p;
          VectorEntry* vse = v.scalar_entries_p;
          if(spv != spuv)  // extend v entry if required
            {
              if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
                {
                  vse=new VectorEntry(*v.scalar_entries_p);
                  std::vector<number_t> dofrenumber;
                  if(u_p == v.u_p) dofrenumber=renumber(v.cdofs(), cdofs());
                  else              dofrenumber=renumber(v.cdofs(), dualDofComponents(cdofs()));
                  if(dofrenumber.size() != 0) vse->extendEntries(dofrenumber, v.cdofs().size());
                }
              if(v.entries_p != 0)
                {
                  ve = new VectorEntry(*v.entries_p);      //copy of v entries
                  std::vector<number_t> dofrenumber = renumber(spuv, spv);
                  if(dofrenumber.size() != 0) ve->extendEntries(dofrenumber, spuv->dimSpace());      //do algebraic operation on entries
                }
            }
          if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
            {
              *scalar_entries_p += *vse;
            }
          if(v.entries_p != 0)
            {
              *entries_p += *ve;  // do algebraic operation on entries
            }
        }
    }

  return *this;
}

SuTermVector& SuTermVector::operator-=(const SuTermVector& v)           // operation U-=V
{
  if(u_p != v.u_p && u_p != v.u_p->dual_p())
    {
      where("SutermVector::operator-=");
      error("term_inconsistent_unknowns"); //different unknowns, error
    }
  if(v.entries_p == 0 && v.scalar_entries_p == 0) return *this;    //v has no entries, add 0 !

  // first, we have to check if both vectors have the same entries structure
  // if no, we have to set v to the same structure, through a copy
  bool copyNeeded=false;
  bool vToScalar=false;
  bool vToVector=false;
  if(scalar_entries_p != 0 && scalar_entries_p != entries_p && v.scalar_entries_p == 0)
    {
      copyNeeded=vToScalar=true;
    }
  else if(scalar_entries_p == 0 && v.scalar_entries_p !=0 && v.scalar_entries_p != v.entries_p)
    {
      copyNeeded=vToVector=true;
    }

  if(copyNeeded)
    {
      SuTermVector vcopy(v);
      if(vToScalar)
        {
          vcopy.toScalar();
        }
      if(vToVector)
        {
          vcopy.toVector();
        }

      Space* spu = space_p, *spv = v.space_p;
      if(spu == spv)  // same numbering => go to entries operation
        {
          if(vToVector)
            {
              if(entries_p == 0) entries_p = new VectorEntry(vcopy.valueType(), vcopy.strucType(), vcopy.nbDofs(), vcopy.nbOfComponents());
              *entries_p -= *vcopy.entries_p;
            }
          if(vToScalar)
            {
              if(scalar_entries_p == 0) scalar_entries_p = new VectorEntry(vcopy.valueType(), vcopy.strucType(), vcopy.nbDofs()*vcopy.nbOfComponents(), 1);
              *scalar_entries_p -= *vcopy.scalar_entries_p;
            }
        }
      else // not the same numbering => dynamic renumbering
        {
          Space* spuv = mergeSubspaces(spu,spv); // create union of subspaces, may be spu or spv or anything else
          space_p=spuv;                          // update subspace member

          if(spu != spuv)  // extend current entry if required
            {
              if(vToVector)
                {
                  if(entries_p != 0 && vcopy.entries_p != 0)
                    {
                      std::vector<number_t> dofrenumber = renumber(spuv, spu);
                      if(dofrenumber.size() != 0) entries_p->extendEntries(dofrenumber, spuv->dimSpace());
                    }
                  VectorEntry* ve = vcopy.entries_p;

                  if(spv != spuv)
                    {
                      if(entries_p != 0 && vcopy.entries_p != 0)
                        {
                          ve = new VectorEntry(*vcopy.entries_p); // copy of v entries
                          std::vector<number_t> dofrenumber = renumber(spuv, spv);
                          if(dofrenumber.size() != 0) ve->extendEntries(dofrenumber, spuv->dimSpace());  // do algebraic operation on entries
                        }
                    }
                  if(entries_p != 0 && vcopy.entries_p != 0)
                    {
                      *entries_p -= *ve;  // do algebraic operation on entries
                    }
                }
              if(vToScalar)
                {
                  if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                    {
                      std::vector<number_t> dofrenumber;
                      if(u_p == vcopy.u_p) dofrenumber=renumber(cdofs(), vcopy.cdofs());
                      else                  dofrenumber=renumber(cdofs(), dualDofComponents(vcopy.cdofs()));
                      if(dofrenumber.size() != 0) scalar_entries_p->extendEntries(dofrenumber, cdofs().size());
                    }
                  VectorEntry* vse = vcopy.scalar_entries_p;

                  if(spv != spuv)
                    {
                      if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                        {
                          vse=new VectorEntry(*vcopy.scalar_entries_p);
                          std::vector<number_t> dofrenumber;
                          if(u_p == vcopy.u_p) dofrenumber=renumber(vcopy.cdofs(), cdofs());
                          else                  dofrenumber=renumber(vcopy.cdofs(), dualDofComponents(cdofs()));
                          if(dofrenumber.size() != 0) vse->extendEntries(dofrenumber, vcopy.cdofs().size());
                        }
                    }
                  if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
                    {
                      *scalar_entries_p -= *vse;  // do algebraic operation on entries
                    }
                }
            }
        }
    }
  else
    {
      Space* spu = space_p, *spv = v.space_p;
      if(spu == spv)  // same numbering => go to entries operation
        {
          if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
            {
              if(scalar_entries_p == 0) scalar_entries_p = new VectorEntry(v.valueType(), v.strucType(), v.nbDofs()*v.nbOfComponents(), 1);
              if(scalar_entries_p != entries_p) *scalar_entries_p -= *v.scalar_entries_p;
            }
          if(v.entries_p != 0)
            {
              if(entries_p == 0) entries_p = new VectorEntry(v.valueType(), v.strucType(), v.nbDofs(), v.nbOfComponents());
              *entries_p -= *v.entries_p;
            }
        }
      else // not the same numbering => dynamic renumbering
        {
          Space* spuv = mergeSubspaces(spu, spv); // create union of subspaces, may be spu or spv or anything else
          space_p = spuv;                         // update subspace member

          if(spu != spuv)  //extend current entry if required
            {
              if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
                {
                  std::vector<number_t> dofrenumber;
                  if(u_p == v.u_p) dofrenumber=renumber(cdofs(), v.cdofs());
                  else              dofrenumber=renumber(cdofs(), dualDofComponents(v.cdofs()));
                  if(dofrenumber.size() != 0) scalar_entries_p->extendEntries(dofrenumber, cdofs().size());
                }
              if(v.entries_p != 0)
                {
                  std::vector<number_t> dofrenumber = renumber(spuv, spu);
                  if(dofrenumber.size() != 0) entries_p->extendEntries(dofrenumber, spuv->dimSpace());
                }
            }
          VectorEntry* ve = v.entries_p;
          VectorEntry* vse = v.scalar_entries_p;
          if(spv != spuv)  // extend v entry if required
            {
              if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
                {
                  vse=new VectorEntry(*v.scalar_entries_p);
                  std::vector<number_t> dofrenumber;
                  if(u_p == v.u_p) dofrenumber=renumber(v.cdofs(), cdofs());
                  else              dofrenumber=renumber(v.cdofs(), dualDofComponents(cdofs()));
                  if(dofrenumber.size() != 0) vse->extendEntries(dofrenumber, v.cdofs().size());
                }
              if(v.entries_p != 0)
                {
                  ve = new VectorEntry(*v.entries_p);      //copy of v entries
                  std::vector<number_t> dofrenumber = renumber(spuv, spv);
                  if(dofrenumber.size() != 0) ve->extendEntries(dofrenumber, spuv->dimSpace());      //do algebraic operation on entries
                }
            }
          if(v.scalar_entries_p != 0 && v.scalar_entries_p != v.entries_p)
            {
              *scalar_entries_p -= *vse;
            }
          if(v.entries_p != 0)
            {
              *entries_p -= *ve;  // do algebraic operation on entries
            }
        }
    }

  return *this;
}

// merging a SuTermVector with current one BUT preserving the current value when common dofs
// like a += but without summing values of common dofs
// process works as following :
//    - detect common dofs
//    - if there are some
//         - copy of SuTermVector to merge with
//         - replace by 0 the value related to a common dofs
//    - use += operator
SuTermVector& SuTermVector::merge(const SuTermVector& v)
{
  trace_p->push("SuTermVector::merge(SuTermVector)");
  if(u_p != v.u_p && u_p != v.u_p->dual_p())
    {
      where("SutermVector::merge");
      error("term_inconsistent_unknowns"); //different unknowns, error
    }
  Space* spu = space_p, *spv = v.space_p;
  if(spu == spv)    //same numbering, merge do nothing
    {
      trace_p->pop();
      return *this;
    }
  //find common dofs
  std::vector<number_t> dofrenumber = renumber(spu, spv);
  bool noCommonDofs = true;
  std::vector<number_t>::iterator itn = dofrenumber.begin();
  while(itn != dofrenumber.end() && noCommonDofs)
    {
      if(*itn != 0) noCommonDofs = false;
      else ++itn;
    }
  SuTermVector* vc = const_cast<SuTermVector*>(&v);
  if(!noCommonDofs)   //have common dofs, copy v and set 0 values on common dofs
    {
      vc = new SuTermVector(v);
      StrucType st = vc->strucType();
      ValueType vt = vc->valueType();
      VectorEntry* entries = vc->entries();
      complex_t z0(0., 0.);
      switch(st)
        {
          case _scalar :
            if(vt == _real)
              {
                number_t i = 1;
                for(itn = dofrenumber.begin(); itn != dofrenumber.end(); ++itn, ++i)
                  if(*itn != 0)  entries->setValue(i, 0.);
              }
            else
              {
                number_t i = 1;
                for(itn = dofrenumber.begin(); itn != dofrenumber.end(); ++itn, ++i)
                  if(*itn != 0)  entries->setValue(i, z0);
              }
            break;
          case _vector :
            if(vt == _real)
              {
                number_t i = 1;
                Vector<real_t> rv0(vc->nbOfComponents(), 0.);
                for(itn = dofrenumber.begin(); itn != dofrenumber.end(); ++itn, ++i)
                  if(*itn != 0)  entries->setValue(i, rv0);
              }
            else
              {
                number_t i = 1;
                Vector<complex_t> rc0(vc->nbOfComponents(), z0);
                for(itn = dofrenumber.begin(); itn != dofrenumber.end(); ++itn, ++i)
                  if(*itn != 0)  entries->setValue(i, rc0);
              }
            break;
          default :
            error("scalar_or_vector");
        }
    }
  //add to current
  (*this) += *vc;

  if(vc != &v) delete vc;
  trace_p->pop();
  return *this;
}

/*! compute SuTermVector from a linear combination of SuTermVector's
 the current SuTermVector is assumed to be void, this function fills every things
 if it is not the case, every thing is clear
*/
void SuTermVector::compute(const LcTerm<SuTermVector>& lc)
{
  trace_p->push("SuTermVector::compute(LcTerm)");
  clear();
  if(lc.size() == 0) error("is_void", "LcTerm");

  termType_ = _sutermVector;
  if(lc.size() == 1)  //copy and multiplication
    {
      const SuTermVector* sut = dynamic_cast<const SuTermVector*>(lc.begin()->first);
      if(sut == 0) error("null_pointer", "SuTermVector");
      complex_t a = lc.begin()->second;
      sulf_p = sut->sulf_p;
      u_p = sut->u_p;
      space_p = sut->space_p;
      subspaces = sut->subspaces;
      entries_p = new VectorEntry(*sut->entries_p);
      if(a != complex_t(1., 0.)) *entries_p *= a;
      computed() = true;
      trace_p->pop();
      return;
    }

  //clear current linear form properties
  subspaces.clear();
  if(sulf_p == 0) delete sulf_p;
  sulf_p = 0;

  // create largest space and subspace numberings related to this largest space
  std::vector<Space*> spaces(lc.size(), 0);
  std::vector<Space*>::iterator its = spaces.begin();
  const SuTermVector* stv = dynamic_cast<const SuTermVector*>(lc.begin()->first);
  if(stv == 0) error("null_pointer", "SuTermVector");
  ValueType vt = stv->valueType();
  StrucType st = stv->strucType();
  u_p = stv->u_p;
  for(LcTerm<SuTermVector>::const_iterator itt = lc.begin(); itt != lc.end(); itt++, its++)
    {
      const SuTermVector* stv = dynamic_cast<const SuTermVector*>(itt->first);
      if(stv == 0) error("null_pointer", "SuTermVector");
      *its = stv->space_p; // list of (sub)spaces
      if(stv->valueType() == _complex) vt = _complex;
      if(stv->strucType() != st)
        error("term_mismatch_structures", words("structure", stv->strucType()), words("structure", st));
    }
  if(lc.coefType() == _complex) vt = _complex;
  space_p = mergeSubspaces(spaces, true);

  //allocate entries of current SuTermVector
  entries_p = new VectorEntry(vt, st, space_p->dimSpace(), u_p->nbOfComponents());

  //do linear combination
  its = spaces.begin();
  for(LcTerm<SuTermVector>::const_iterator itt = lc.begin(); itt != lc.end(); itt++, its++)
    {
      const SuTermVector* stv = dynamic_cast<const SuTermVector*>(itt->first);
      if(!stv->computed()) error("not_computed_term", stv->name());
      Space* sp = *its;
      //if (sp->typeOfSpace() == _subSpace)
      if(sp->typeOfSpace() == _subSpace && sp != space_p)
        {
          if(vt == _real) entries_p->add(*stv->entries_p, sp->subSpace()->dofNumbers(), itt->second.real());
          else entries_p->add(*stv->entries_p, sp->subSpace()->dofNumbers(), itt->second);
        }
      else  //subspace is the largest space (same numbering)
        {
          if(vt == _real) entries_p->add(*stv->entries_p, itt->second.real());
          else entries_p->add(*stv->entries_p, itt->second);
        }
    }
  //may be it is best to delete new subspaces created here (to be done later)
  computed() = true;
  trace_p->pop();
}

/*! interpolate current sutermvector on u space and domain dom
           sutr(i) = sum_j sutv_j w_j(Mi)  Mi nodes of Lagrange u-space, Mi in dom
                                           wj shape functions related to sutv
    return a new SuTermMatrix allocated here
    quite expansive !!!
*/
SuTermVector* SuTermVector::interpolate(const Unknown& u, const GeomDomain& dom)
{
  trace_p->push("SuTermVector::interpolate(Unknown, GeomDomain)");
  if(u.space() == space_p) return onDomain(dom);  //same space, apply restriction to domain dom
  if(!u.space()->isFE()) error("not_fe_space_type", u.space()->name());
  if(u.space()->interpolation()->type != _Lagrange) error("lagrange_fe_space_only", u.space()->name());

  //get subspace
  Space* spu = u.space();
  Space* subsp = Space::findSubSpace(&dom, spu);         //find subspace (linked to domain dom) of spu
  if(subsp == 0) subsp = new Space(dom, *spu, spu->name() + "_" + dom.name());  //create if not exist

  //create suTermVector
  std::vector<number_t> dofs = subsp->dofIds();
  number_t nd = dofs.size();
  string_t na = name_ + "_interpolated_" + u.space()->name() + "_" + dom.name();
  SuTermVector* sutv = new SuTermVector(na, &u, subsp, valueType(), nd, u.nbOfComponents());

  //travel dofs and interpolate at node, store value
  std::vector<number_t>::iterator itd = dofs.begin();
  number_t k = 1;
  for(; itd != dofs.end(); ++itd, ++k)
    {
      Point p = subsp->rootSpace()->feSpace()->dofs[*itd - 1].coords();
      Value v = evaluate(p);
      sutv->entries_p->setValue(k, v);
    }
  sutv->computed() = true;
  trace_p->pop();
  return sutv;

}

//! print SuTermVector
void SuTermVector::print(std::ostream& out) const
{
  print(out, false, true);
}

//if raw=true only x y z v is print
void SuTermVector::print(std::ostream& out, bool raw, bool header) const
{
  if(theVerboseLevel == 0) return;
  if(header && !raw) out << "   SuTermVector " << name() << ", " << " unknown " << u_p->name() << ", "
                           << words("value", valueType()) << " " << words("structure", strucType())
                           << ", ";
  if(sulf_p != 0) out << (*sulf_p) << "\n";
  if(!computed())
    {
      out << words("not computed") << "\n";
      return;
    }
  if(entries_p == 0 && scalar_entries_p == 0)
    {
      out << "  unallocated SuTermVector block" << eol;
      return;
    }
  if(theVerboseLevel < 2) return;
  if(!raw)
    {
      if(nbOfComponents() > 1) out << " (size=" << nbOfComponents() << "),";
      out << " " << size() << " " << words("values");
      if(theVerboseLevel < 2)
        {
          out << eol;
          return;
        }

      //first and last values
      blanks(eol, 3);
      number_t s = size();
      number_t n = std::min(number_t(theVerboseLevel), s);

      if(entries_p != 0)  //display dofid and coordinates when dof is a ponctual fedof
        {
          number_t dims = space_p->nbDofs();
          number_t wd = std::floor(std::log10(dims))+1;
          if(2 * n >= dims)  //print all value
            {
              out << ", all values" << eol;
              std::vector<string_t> com(dims, "");
              number_t w=0;
              for(number_t i = 0; i < dims; i++)
                { com[i] = space_p->dof(i + 1).dofData(wd); w=std::max(w,com[i].size());}
              for(number_t i = 0; i <dims; i++)
                {
                  number_t s = w-com[i].size();
                  if(s>0) com[i]+= string_t(s,' ');
                  com[i]+=" -> ";
                }
              entries_p->printFirst(out, dims, com);
            }
          else
            {
              out << ", first and last " << tostring(n) << " values" << eol;
              std::vector<string_t> comd(n, ""), comf(n, "");
              number_t w=0;
              for(number_t i = 0; i < n; i++)
                {
                  comd[i] = space_p->dof(i + 1).dofData(wd); w=std::max(w,comd[i].size());
                  comf[i] = space_p->dof(dims - n + i + 1).dofData(wd); w=std::max(w,comf[i].size());
                }
              for(number_t i = 0; i <n; i++)
                {
                  number_t s = w-comd[i].size();
                  if(s>0) comd[i]+= string_t(s,' ');
                  comd[i]+=" -> ";
                  s = w-comf[i].size();
                  if(s>0) comf[i]+= string_t(s,' ');
                  comf[i]+=" -> ";
                }
              entries_p->printFirst(out, n, comd);
              out << "..." << eol;
              entries_p->printLast(out, n, comf);
            }
        }
      if(scalar_entries_p != 0)
        {
          out << " scalar entries : " << tostring(scalar_entries_p->size()) << " "
              << words("value", valueType()) << " values" << eol;
          out << words("firste") << " " << words("row") << " cdofs : ";
          number_t nbd = cdofs_.size();
          number_t m = std::min(number_t(theVerboseLevel), nbd);
          std::vector<DofComponent>::const_iterator itd = cdofs_.begin();
          for(number_t i = 0; i < m; i++, itd++)  out  << *itd << "  ";
          if(m < nbd) out << "...";
          out << eol;
          if(scalar_entries_p != entries_p) out << (*scalar_entries_p);
        }
      blanks(eol, -3);
      out << eol;
    }
  else  //raw format
    {
      if(entries_p != 0)
        {
          number_t dims = space_p->nbDofs();
          std::vector<string_t> com(dims, "");
          for(number_t i = 0; i < dims; i++)
            {
              Point pi = space_p->dof(i + 1).coords();
              for(number_t j = 0; j < pi.size(); j++) com[i] +=  tostring(pi[j]) + " ";
            }
          entries_p->printFirst(out, dims, com);
        }
      if(scalar_entries_p != 0 && scalar_entries_p != entries_p)
        {
          number_t dims = nbOfComponents() * space_p->nbDofs();
          std::vector<string_t> com(dims, "");
          scalar_entries_p->printFirst(out, dims, com);
        }
    }
}

complex_t innerProduct(const SuTermVector& tv1, const SuTermVector& tv2)
{
  trace_p->push("innerProduct(SuTermVector,SuTermVector)");
  if(&tv1 == &tv2)  //same SuTermVectors
    {
      if(tv1.scalar_entries() != 0)
        {
          trace_p->pop();
          return innerProduct(*tv1.scalar_entries(), *tv1.scalar_entries());
        }
      else
        {
          trace_p->pop();
          return innerProduct(*tv1.entries(), *tv1.entries());
        }
    }
  //check consistancy of unknowns
  complex_t res = 0.;
  if(tv1.up() != tv2.up() && tv1.up() != tv2.up()->dual_p()) error("term_inconsistent_unknowns");
  number_t n1 = tv1.spacep()->dimSpace(), n2 = tv2.spacep()->dimSpace();
  if(n1 > n2) res = innerProduct(tv2, tv1);
  else
    {
      bool copyNeeded2=false;
      bool toScalar2=false;
      bool toVector2=false;
      if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() && tv2.scalar_entries() == 0)
        {
          copyNeeded2=toScalar2=true;
        }
      else if(tv1.scalar_entries() == 0 && tv2.scalar_entries() !=0 && tv2.scalar_entries() != tv2.entries())
        {
          copyNeeded2=toVector2=true;
        }

      if(copyNeeded2)
        {
          SuTermVector tv2copy(tv2);
          if(toScalar2) tv2copy.toScalar();
          if(toVector2) tv2copy.toVector();
          if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() &&
              tv2copy.scalar_entries() != 0 && tv2copy.scalar_entries() != tv2copy.entries())
            {
              if(tv1.spacep() == tv2copy.spacep()) res = innerProduct(*tv1.scalar_entries(), *tv2copy.scalar_entries());
              else  //align vector
                {
                  bool useDual = (tv1.up() == tv2.up()->dual_p());
                  tv2copy.extendScalarTo(tv1.cdofs(), useDual);      //align tv2copy on tv1
                  res = innerProduct(*tv1.scalar_entries(), *tv2copy.scalar_entries());
                }
            }
          else
            {
              if(tv1.entries() == 0) error("null_pointer", "tv1.entries_p");
              if(tv2copy.entries() == 0) error("null_pointer", "tv2.entries_p");
              if(tv1.spacep() == tv2copy.spacep()) res = innerProduct(*tv1.entries(), *tv2copy.entries());
              else  //align vector
                {
                  tv2copy.extendTo(tv1);      //align tv2copy on tv1
                  res = innerProduct(*tv1.entries(), *tv2copy.entries());
                }
            }
        }
      else
        {
          if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() &&
              tv2.scalar_entries() != 0 && tv2.scalar_entries() != tv2.entries())
            {
              if(tv1.spacep() == tv2.spacep()) res = innerProduct(*tv1.scalar_entries(), *tv2.scalar_entries());
              else  //align vector
                {
                  SuTermVector tw2(tv2);  //copy tv2
                  bool useDual = (tv1.up() == tv2.up()->dual_p());
                  tw2.extendScalarTo(tv1.cdofs(), useDual);      //align tw2 on tv1
                  res = innerProduct(*tv1.scalar_entries(), *tw2.scalar_entries());
                }
            }
          else
            {
              if(tv1.entries() == 0) error("null_pointer", "tv1.entries_p");
              if(tv2.entries() == 0) error("null_pointer", "tv2.entries_p");
              if(tv1.spacep() == tv2.spacep()) res = innerProduct(*tv1.entries(), *tv2.entries());
              else  //align vector
                {
                  SuTermVector tw2(tv2);  //copy tv2
                  tw2.extendTo(tv1);      //align tw2 on tv1
                  res = innerProduct(*tv1.entries(), *tw2.entries());
                }
            }
        }
    }
  trace_p->pop();
  return res;
}

complex_t hermitianProduct(const SuTermVector& tv1, const SuTermVector& tv2)
{
  trace_p->push("hermitianProduct(SuTermVector,SuTermVector)");
  complex_t res = 0;
  if(&tv1 == &tv2)  //same SuTermVectors
    {
      if(tv1.scalar_entries() != 0)
        {
          trace_p->pop();
          return hermitianProduct(*tv1.scalar_entries(), *tv1.scalar_entries());
        }
      else
        {
          trace_p->pop();
          return hermitianProduct(*tv1.entries(), *tv1.entries());
        }
    }
  //check consistancy of unknowns
  if(tv1.up() != tv2.up() && tv1.up() != tv2.up()->dual_p()) error("term_inconsistent_unknowns");
  number_t n1 = tv1.spacep()->dimSpace(), n2 = tv2.spacep()->dimSpace();
  if(n1 > n2) res = std::conj(hermitianProduct(tv2, tv1));
  else
    {
      bool copyNeeded2=false;
      bool toScalar2=false;
      bool toVector2=false;
      if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() && tv2.scalar_entries() == 0)
        {
          copyNeeded2=toScalar2=true;
        }
      else if(tv1.scalar_entries() == 0 && tv2.scalar_entries() != 0 && tv2.scalar_entries() != tv2.entries())
        {
          copyNeeded2=toVector2=true;
        }

      if(copyNeeded2)
        {
          SuTermVector tv2copy(tv2);
          if(toScalar2) tv2copy.toScalar();
          if(toVector2) tv2copy.toVector();
          if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() &&
              tv2copy.scalar_entries() != 0 && tv2copy.scalar_entries() != tv2copy.entries())
            {
              if(tv1.spacep() == tv2copy.spacep()) res = hermitianProduct(*tv1.scalar_entries(), *tv2copy.scalar_entries());
              else  //align vector
                {
                  bool useDual = (tv1.up() == tv2.up()->dual_p());
                  tv2copy.extendScalarTo(tv1.cdofs(), useDual);      //align tv2copy on tv1
                  res = hermitianProduct(*tv1.scalar_entries(), *tv2copy.scalar_entries());
                }
            }
          else
            {
              if(tv1.entries() == 0) error("null_pointer", "tv1.entries_p");
              if(tv2copy.entries() == 0) error("null_pointer", "tv2.entries_p");
              if(tv1.spacep() == tv2copy.spacep()) res = hermitianProduct(*tv1.entries(), *tv2copy.entries());
              else  //align vector
                {
                  tv2copy.extendTo(tv1);      //align tv2copy on tv1
                  res = hermitianProduct(*tv1.entries(), *tv2copy.entries());
                }
            }
        }
      else
        {
          if(tv1.scalar_entries() != 0 && tv1.scalar_entries() != tv1.entries() &&
              tv2.scalar_entries() != 0 && tv2.scalar_entries() != tv2.entries())
            {
              if(tv1.spacep() == tv2.spacep()) res = hermitianProduct(*tv1.scalar_entries(), *tv2.scalar_entries());
              else  //align vector
                {
                  SuTermVector tw2(tv2);  //copy tv2
                  bool useDual = (tv1.up() == tv2.up()->dual_p());
                  tw2.extendScalarTo(tv1.cdofs(), useDual);      //align tw2 on tv1
                  res = hermitianProduct(*tv1.scalar_entries(), *tw2.scalar_entries());
                }
            }
          else
            {
              if(tv1.entries() == 0) error("null_pointer", "tv1.entries_p");
              if(tv2.entries() == 0) error("null_pointer", "tv2.entries_p");
              if(tv1.spacep() == tv2.spacep()) res = hermitianProduct(*tv1.entries(), *tv2.entries());
              else  //align vector
                {
                  SuTermVector tw2(tv2);  //copy tv2
                  tw2.extendTo(tv1);      //align tw2 on tv1
                  res = hermitianProduct(*tv1.entries(), *tw2.entries());
                }
            }
        }
    }
  trace_p->pop();
  return res;
}

// vector norms
real_t norm1(const SuTermVector& tv)
{
  if(tv.scalar_entries() != 0) return norm1(*tv.scalar_entries());
  if(tv.entries() != 0) return norm1(*tv.entries());
  return 0.;
}

real_t norminfty(const SuTermVector& tv)
{
  if(tv.scalar_entries() != 0) return norminfty(*tv.scalar_entries());
  if(tv.entries() != 0) return norminfty(*tv.entries());
  return 0.;
}

// return the value (in complex) of component being the largest one in absolute value
// travel all SCALAR components of SuTermVector using scalar entries if exists
complex_t SuTermVector::maxValAbs() const
{
  if(scalar_entries_p != 0) return scalar_entries_p->maxValAbs();
  if(entries_p != 0) return entries_p->maxValAbs();
  return complex_t(0.);
}

//=================================================================================
// SuTermVectors member functions
//=================================================================================

// move all SuTermVector's to real
SuTermVectors& SuTermVectors::toReal()
{
  std::vector<SuTermVector>::iterator itv = begin();
  for(; itv != end(); itv++) itv->toReal();
  return *this;
}

SuTermVectors& SuTermVectors::toImag()
{
  std::vector<SuTermVector>::iterator itv = begin();
  for(; itv != end(); itv++) itv->toImag();
  return *this;
}

void SuTermVectors::print(std::ostream& os) const
{
  os << "list of SuTermVector's (" << size() << ")" << eol;
  if(theVerboseLevel < 2) return;
  std::vector<SuTermVector>::const_iterator itv = begin();
  for(; itv != end(); itv++) itv->print(os);
}

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

/*! adjust the  scalar entries to the cdofs numbering newcdofs
    when cdofs  is included in newcdofs, new zeros are introduced (extension)
    when newcdofs is included in cdofs , coefficients are removed (restriction)
    else it is in a same time a restriction (removed cdofs)  and an extension (omitted cdof)
    in any case, the coefficient may be permuted and cdofs = newcdofs at the end

    used by SuTermVector and TermVector
*/
void adjustScalarEntriesG(VectorEntry*& scalar_entries_p, std::vector<DofComponent>& cdofs, const std::vector<DofComponent>& newcdofs)
{
  trace_p->push("adjustScalarEntriesG(VectorEntry*, vector<DofComponent>&, vector<DofComponent>&)");
  std::map<DofComponent, number_t> mcdofs;
  std::vector<DofComponent>::iterator itd = cdofs.begin();
  std::vector<DofComponent>::const_iterator itn = newcdofs.begin();
  number_t k = 0;
  bool sameNum = true;
  for(; itd != cdofs.end(); itd++, k++)
    {
      mcdofs[*itd] = k;
      if(sameNum && itn != newcdofs.end())
        {
          if(*itd != *itn  && *itd != itn->dual()) sameNum = false;  //allows dual unknown (permissive)
          itn++;
        }
    }
  if(sameNum && cdofs.size() == newcdofs.size())
    {
      trace_p->pop();  //nothing done
      return;
    }

  //create new vector
  ValueType vt = _real;
  if(scalar_entries_p != 0) vt = scalar_entries_p->valueType_;
  VectorEntry* ve = new VectorEntry(vt, _scalar, newcdofs.size());
  std::map<DofComponent, number_t>::iterator itmd;
  if(vt == _real)  //real case
    {
      Vector<real_t>::iterator itve = ve->beginr(), itv = scalar_entries_p->beginr();
      for(itn = newcdofs.begin(); itn != newcdofs.end(); itn++, itve++)
        {
          itmd = mcdofs.find(*itn);
          if(itmd == mcdofs.end()) itmd = mcdofs.find(itn->dual());
          if(itmd != mcdofs.end()) *itve = *(itv + itmd->second);
        }
    }
  else  //complex case
    {
      Vector<complex_t>::iterator itve = ve->beginc(), itv = scalar_entries_p->beginc();
      for(itn = newcdofs.begin(); itn != newcdofs.end(); itn++, itve++)
        {
          itmd = mcdofs.find(*itn);
          if(itmd == mcdofs.end()) itmd = mcdofs.find(itn->dual());  //find with dual (permissive)
          if(itmd != mcdofs.end()) *itve = *(itv + itmd->second);
        }
    }
  //update current vector
  if(scalar_entries_p != 0) delete scalar_entries_p;
  scalar_entries_p = ve;
  cdofs = newcdofs;
  trace_p->pop();
}

} //end of namespace xlifepp